弹道模型matlab代码

标题：探索MATLAB中各类拟合曲线的代码应用在MATLAB中，拟合曲线是数据分析和模型建立中常用的技术之一。

通过拟合曲线，我们可以了解数据之间的关联性并建立预测模型，为进一步分析和应用数据奠定基础。

本文将深入探讨MATLAB中各类拟合曲线的代码应用，帮助读者更深入地理解该主题。

一、线性拟合曲线1. 使用MATLAB进行线性拟合曲线的代码示例在MATLAB中，使用polyfit函数可以进行线性拟合。

对一组数据点(x, y)进行线性拟合，代码如下：```matlabx = [1, 2, 3, 4, 5];y = [2, 3.5, 5, 7, 8.5];p = polyfit(x, y, 1);```其中，x为自变量，y为因变量，1表示进行一次线性拟合。

通过polyfit函数，可以得到线性拟合的系数p。

2. 线性拟合曲线的应用和特点线性拟合曲线适用于线性关系较为明显的数据，例如物理实验数据中的直线关系。

通过线性拟合，可以获得各项系数，对数据进行预测和建模。

二、多项式拟合曲线1. 使用MATLAB进行多项式拟合曲线的代码示例在MATLAB中，使用polyfit函数同样可以进行多项式拟合。

对一组数据点(x, y)进行二次多项式拟合，代码如下：```matlabx = [1, 2, 3, 4, 5];y = [1, 4, 9, 16, 25];p = polyfit(x, y, 2);```其中，x为自变量，y为因变量，2表示进行二次多项式拟合。

通过polyfit函数，同样可以得到多项式拟合的系数p。

2. 多项式拟合曲线的应用和特点多项式拟合曲线适用于数据中存在曲线关系的情况，通过选择合适的最高次数，可以灵活地拟合各种曲线形状。

三、非线性拟合曲线1. 使用MATLAB进行非线性拟合曲线的代码示例在MATLAB中，使用fit函数可以进行非线性拟合。

对一组数据点(x, y)进行指数函数拟合，代码如下：```matlabx = [1, 2, 3, 4, 5];y = [2.1, 7.4, 16.1, 29.3, 48.2];f = fit(x', y', 'exp1');```其中，x为自变量，y为因变量，'exp1'表示进行指数函数拟合。

解弹道方程求解舰炮武器系统射击诸元的数学模型一、引言随着现代军事技术的发展，舰炮武器系统作为海军舰船上重要的作战装备之一，其射击精度和命中率直接影响了作战效果。

在舰炮武器系统的设计和研发过程中，需要进行射击参数的分析和优化，而解弹道方程是求解舰炮武器系统射击诸元的数学模型中的重要内容之一。

二、解弹道方程的基本原理解弹道方程是舰炮武器系统射击诸元的数学模型中的关键环节，其基本原理是根据枪弹在发射过程中受到的空气阻力、重力等外力的影响，建立弹道方程并求解得到枪炮射程、射速和射击精度等参数。

三、弹道方程的建立1. 枪弹运动方程枪弹在发射过程中受到重力、空气阻力等外力的作用，运动方程可以描述为：弹道方程1弹道方程2弹道方程32. 外力分析在建立弹道方程的过程中，需要对枪弹受到的外力进行分析，例如重力、空气阻力等的影响。

3. 初值条件在求解弹道方程时，需要确定一定的初值条件，包括发射角度、初速度以及大气条件等。

四、弹道方程的求解1. 数值计算方法弹道方程通常采用数值方法进行求解，例如有限差分法、龙格-库塔法等数值计算方法。

2. 计算工具为了求解弹道方程，通常需要借助专业的计算工具，例如MATLAB、Python等。

3. 求解结果分析求解得到弹道方程后，需要对求解结果进行分析，包括射程、射速、命中率等参数的评估。

五、应用与优化1. 已有案例分析结合实际的舰炮武器系统，可以进行已有案例的弹道方程求解及优化分析，以提高射击精度。

2. 参数优化通过对弹道方程中的各项参数进行优化，可以提高舰炮武器系统的射击精度和命中率。

六、结论解弹道方程求解舰炮武器系统射击诸元的数学模型，对于提高舰炮武器系统的射击精度和命中率具有重要意义。

通过建立弹道方程、求解和优化，能够有效提升舰炮武器系统的作战效能，为我国海军现代化建设做出重要贡献。

七、参考文献1. 《舰炮武器系统原理与设计》，XXX，XXX出版社，2000.2. 《武器装备使用手册》，XXX，XXX出版社，2010.以上就是关于解弹道方程求解舰炮武器系统射击诸元的数学模型的文章内容，希望对读者有所帮助。

1.绘制云图Ex=18En=2He=0.2hold onfor i=1:1000Enn=randn(1)*He+En;x(i)=randn(1)*Enn+Ex;y(i)=exp(-(x(i)-Ex)^2/(2*Enn^2)); plot(x(i),y(i),'*')endEx=48.7En=9.1He=0.39hold onfor i=1:1000Enn=randn(1)*He+En;x(i)=randn(1)*Enn+Ex;y(i)=exp(-(x(i)-Ex)^2/(2*Enn^2)); plot(x(i),y(i),'*')end2.求期望、熵及超熵X1=[51.93 52.51 54.70 43.14 43.85 44.48 44.61 52.08];Y1=[0.91169241573 0.921875 0.96032303371 0.75737359551 0.76983848315 0.7808988764 0.78318117978 0.9143258427];m=8;Ex=mean(X1)En1=zeros(1,m);for i=1:mEn1(1,i)=abs(X1(1,i)-Ex)/sqrt(-2*log(Y1(1,i)));endEn=mean(En1);He=0;for i=1:mHe=He+(En1(1,i)-En)^2;endEn=mean(En1)He=sqrt(He/(m-1))3.平顶山so2环境:X1=[0.013 0.04 0.054 0.065 0.07 0.067 0.058 0.055 0.045];Y1=[0.175675676 0.540540541 0.72972973 0.8783783780.945945946 0.905405405 0.783783784 0.743243243 0.608108108]; m=9;Ex=mean(X1)En1=zeros(1,m);for i=1:mEn1(1,i)=abs(X1(1,i)-Ex)/sqrt(-2*log(Y1(1,i)));endEn=mean(En1);He=0;for i=1:mHe=He+(En1(1,i)-En)^2;endEn=mean(En1)He=sqrt(He/(m-1))1.绘制正向云图Ex=18En=2He=0.2hold onfor i=1:1000Enn=randn(1)*He+En;x(i)=randn(1)*Enn+Ex;y(i)=exp(-(x(i)-Ex)^2/(2*Enn^2));plot(x(i),y(i),'*')endEx=48.7En=9.1He=0.39hold onfor i=1:1000Enn=randn(1)*He+En;x(i)=randn(1)*Enn+Ex;y(i)=exp(-(x(i)-Ex)^2/(2*Enn^2));plot(x(i),y(i),'*')end2.逆向云发生器中需要剔除隶属度大于0. 9999 的云滴，剩下个云滴。

DDF模型是一种广泛应用于地震勘探和地球物理探测领域的数学模型，它可以帮助我们更好地理解地下构造和地震波传播规律。

在实际工程中，使用Matlab编程对DDF模型进行仿真和分析是非常常见的。

本文将介绍DDF模型的原理及其在Matlab程序中的实现。

一、DDF模型的原理1. 地震波传播原理在地球物理勘探中，地震波的传播是一项重要的研究内容。

地震波在地下介质中传播时会发生折射、反射和衍射等现象，这些现象受到介质物性的影响。

2. DDF模型概述DDF模型是一种基于弹性波动方程的数学模型，它可以描述地震波在介质中的传播过程。

DDF模型考虑了介质的弹性性质和几何形态，能够较准确地模拟地震波在复杂介质中的传播情况。

3. DDF模型的理论基础DDF模型基于弹性波动方程推导而来，其具体数学表达为一组偏微分方程。

通过对介质物性和边界条件的合理假设，可以得到DDF模型的数值解，从而实现对地震波传播的模拟和分析。

二、DDF模型的Matlab编程实现1. 编程环境准备在进行DDF模型的Matlab编程之前，首先需要准备好编程环境。

包括安装Matlab软件、了解Matlab的基本语法和数据处理方法等。

2. DDF模型的数值求解DDF模型的数值求解是整个Matlab编程过程的核心部分。

通过将DDF模型的偏微分方程离散化，可以得到一个关于介质物性和边界条件的代数方程组，利用Matlab的数值计算能力可以求解这组方程。

3. 结果可视化在得到DDF模型的数值解之后，还需要对模拟结果进行可视化处理。

可以利用Matlab的绘图功能，将地震波的传播情况以图像的形式清晰展现出来，便于工程人员进行分析和理解。

三、DDF模型在地震勘探中的应用实例1. 地震成像DDF模型在地震成像中有着广泛的应用。

通过对地震波在不同介质中的传播情况进行模拟，可以确定地下各层的构造和性质，为地球物理勘探提供重要的参考信息。

2. 地震波的反演地震波反演是地球物理勘探中的重要技术手段，可以通过对地震波在地下介质中的传播进行模拟，反推出地下介质的物性参数。

西工大制导控制系统仿真技术 试验 实验二 导弹的纵向运动仿真建模一、实验目的与要求1．加深对导弹的纵向模型的理解；2．掌握Matlab/Simulink 仿真模型的建模步骤和注意事项； 3．在Matlab/SimuLink 环境下，进行导弹的纵向运动仿真建模；二、实验内容1、首先给出导弹的的简化纵向运动方程：◆ 质心动力学方程：cos sin sin cos dVmP X mg dtd mV P Y mg dt αθθαθ=--=+-◆ 绕质心转动转动动力学方程： zzz d J M dtω= ◆ 质心运动学方程：cos sin dxV dtdyV dt θθ==◆ 绕质心转动动力学方程： z d dtϑω= ◆ 补充几何关系方程： ϑαθ=+◆ 过载计算方程：cos y V d n g dtθθ=+ 式中： m 为导弹当前的质量；V 为导弹速度； P 为导弹推力； X 为导弹受到的阻力； Y 为导弹受到的升力；g 为重力加速度；z J 为导弹的z 轴的转动惯量；z ω为导弹俯仰角速度；y n 为导弹的法向过载；z M 为导弹受到气动俯仰力矩；ϑ为俯仰角；θ为弹道倾角； α为攻角。

2、下面给出导弹的弹体参数： ● 弹体参数：弹体质量：0m =250；转动惯量：315Jz =；参考面积：0.45S =参考长度： 2.5L =；发动机推力：1000P =；质量秒消耗量：0.46dm = ● 阻力计算公式：x X C qS =，其中：220x x x C C C αα=+，00.2x C =，20.0005x C α=● 升力计算公式：y Y C qS =，其中：z y y y z C C C δααδ=+，0.25y C α=，0.05z y C δ=● 俯仰力矩计算公式：z z M m qSL =，其中：z z z z z m m m δααδ=+，-0.1z m α=，-0.024z z m δ=3、下面给出导弹的控制系统：图 1 导弹的俯仰控制系统图中，控制系统的参数：Kwz = -0.05,Kpny = 0.1,Kiny = 2。

基于人大金仓数据库和MATLAB的弹道目标数据管理作者：蔡红军来源：《数字技术与应用》2017年第11期摘要：由于弹道目标的特殊性，真实弹道目标数据显得十分珍贵，如何对历史数据进行保存和管理，便于进行分析，充分发挥试验数据的价值，变得越来越重要。

本文针对这一问题，结合弹道目标数据的特点，提出用人大金仓数据库进行数据管理，同时利用MATLAB软件进行数据分析，即解决了弹道目标数据的存储问题，同时又利用MATLAB软件实现对数据的可视化分析。

关键词：数据库；人大金仓；数据库访问中图分类号：TP311.13 文献标识码：A 文章编号：1007-9416（2017）11-0104-021 引言弹道目标数据由于保密的原因，很难获得真实目标的实测数据，因此弹道目标数据显得十分珍贵。

随着参与外场试验的次数越来越多，逐渐积累了大量宝贵的原始数据，然而缺乏对现有的数据进行梳理，难以发挥有限试验数据的价值，同时随着时间的流逝，关键数据信息难免遗漏。

因此，急需一个数据库对现有的数据进行管理。

MATLAB软件擅长于数据分析和可视化，同时MATLAB软件中有一个为关系型数据库和MATLAB进行数据交互提供接口函数支持的工具包Database Toolbox，其使用简单，功能强大，通过其提供的接口函数极大简化了MATLAB与关系型数据库间的数据交互，而且它支持同时与多个数据库的并行任务和大量数据集的分段输入。

在WINDOWS环境下，Database Toolbox支持JDBC驱动和ODBC驱动，然而Database Toolbox并不支持非结构化大型对象（Large Object，LOB）数据的存储。

本文结合弹道目标数据的特点，采用人大金仓数据库（KingbaseES）和MATLAB软件完成数据的管理和分析，使用KingbaseES中TEXT数据类型，代替LOB数据类型。

本文中的驱动配置方法和代码程序均在MATLAB 2012b版本和KingbaseES V7版本下进行。

航天飞行动力学作业报告——有翼导弹飞行方案和稳定性分析一、问题描述：1.在给定的条件下，计算纵向理想弹道，并给出采用瞬时平衡假设0zz z z m m δααδ+=时所有纵向参数随时间的变化曲线。

2.不考虑气动力下洗影响，以第一问得出的弹道为基础，选取并计算作为特性点的5个以上点处的纵向短周期扰动运动的动力系数，并分析其在特性点处的自由扰动的稳定性，以及计算在各个特性点处弹体传递函数(),(),()y n W s W s W s αδδϑδ 。

二、模型建立：根据给出的飞行条件进行初步分析，可给出如下假设和简化： 1、 近似认为导弹绕弹体轴的转动是无惯性的。

2、 近似认为导弹控制系统理想工作，既无误差，也无时间延迟。

3、 近似认为各种干扰因素对导弹无任何影响。

4、由于侧向运动参数与x 与y 方向舵偏角都是小量，因此可近似认为相关参数可以忽略。

5、近似认为导弹在某个铅锤面内飞行，即其飞行弹道与铅锤面内的弹道差别不大。

6、近似认为俯仰操纵机构的偏转仅取决于纵向运动参数；偏航、滚转操纵机构的偏转仅取决于侧向运动参数。

根据以上假设，我们可以简化得到以下方程组： 质心移动的动力学方程：mmdddddddd =PPPPPPPPαα−XX −mmmmPPmm mm θθ mmdd ddθθdddd =PPPPmm mm αα+YY −mmmmPPPPPPθθ质心移动的运动学方程：dddddddd =Vcos θ dddddddd =Vsin θ 质量方程：ddmmdddd=mm 纵向平衡关系式：0zz z z m m δααδ+=控制方程：14=0=0εε上式适用于全阶段的飞行方案，但是因为每个阶段的参数会有所不同，因此在不同阶段该方程组会有不同的形式，再根据每个阶段的具体的公式进行数值积分就能够得到最终各参数的变化情况。

三、求解弹道1.第一阶段：给定高度导弹释放后，在第一阶段做无动力滑翔，采用给定高度的飞行方案，其控制系统方程有表达式如下：***2000cos(0.000314 1.1)5000(-)+(-)zH x k H H k H H ϕϕδ=×××+=×× 值得注意的是，控制方程中包含开环增益系数，其值的选取关系到在控制系统下的飞行弹道与给定弹道的相合程度，通过matlab 进行循环迭代调试选取使弹道最相合且震荡最微弱的参数k k φφ ，得到阶段一各参数随时间变化的关系如下图所示：图 3-1 第一阶段飞行参数变化曲线图3-1给出了阶段一导弹速度、弹道倾角、导弹质量、水平位移、导弹高度、z 向舵偏角随时间变化的关系，其中各单位分别为m/s 、rad 、kg 、m 、m 、rad ，时间单位为s. 2.第二阶段：等高飞行导弹在水平位移为9100m 时，发动机开始点火，转入水平飞行模式。

第三章 弹道轨迹仿真求解在弹道学中，知道弹道的轨迹是十分重要的，根据弹道轨迹可以求出弹道参数，而求弹道轨迹曲线的过程也就是解弹道微分方程组的过程。

本章将对解弹道微分方程组的求解进行论述。

3.1 理论基础3.1.1 数学工具微分方程是数学科学联系实际问题的主要桥梁之一，它是含有未知函数及其导数的方程。

在工程实际与科学研究中遇到的微分方程往往比较复杂，在很多情况下，都不能给出解析表达式，这些情况下不适宜采用解析法来求解，而需采用数值解法来求近似解[8]。

求解微分方程的问题大体上可以分为初值问题和边值问题两大部分。

关于初值问题的求解，常用的解法有：龙格库塔方法、泰勒级数展开法、外推法、欧拉方程等；边值问题的常用解法有：配置法、有限差分法、打靶法等。

弹道微分方程组的求解属于初值问题的求解，我们已经知道微分方程组的数值解法有很多种，但是在众多的算法中，四阶龙格-库塔法具有比较高的精确度，是在求解微分方程组数值解过程中的一种优先选取的算法。

下面将对四阶龙格-塔法的原理进行简单叙述。

3.1.2 龙格-库塔法原理龙格-库塔（Runge-Kutta ）方法是一种在工程上应用广泛的“高精度单步算法”。

Euler 公式可改写成⎩⎨⎧+==+Ky y y x hf K i i i i 1),( (3-1)则1+i y 的表达式与)(1+i x y 的Taylor 展开式的前两项完全相同，即局部截断误差为)(2h O 。

同理，改进Euler 公式，将其改写成 2112121K K y y i i ++=+ (3-2) 其中),(1i i y x hf K =，),(2h y h x hf K i i ++=。

上述两组公式在形式上共同点：都是用),(y x f 在某些点上值的线性组合得出)(1+i x y 的近似值1+i y ，且增加计算),(y x f 的次数，可提高截断误差的阶。

如欧拉法：每步计算一次),(y x f 的值，为一阶方法。

matlab点的轨迹
在Matlab中，点的轨迹通常由一系列坐标点组成。

每个坐标点都有一个x和y坐标值。

通过在坐标系中连接这些点，就可以绘制出点的轨迹。

例如，可以使用以下代码生成一个包含5个坐标点的轨迹：
matlab复制代码:
x = [1, 2, 3, 4, 5];
y = [1, 4, 9, 16, 25];
plot(x, y);
上述代码将生成一个包含5个坐标点的轨迹，坐标点的x坐标分别为1、2、3、4、5，y坐标分别为1、4、9、16、25。

如果想要绘制动态的点的轨迹，可以使用Matlab中的动态图形工具包（Dynamic Dataflow Blockset）来实现。

这个工具包允许用户通过连接不同的数据流块来创建动态系统，从而模拟点的轨迹。

具体实现方法可以参考Matlab官方文档和教程。

多元midas模型和单变量midas的matlab代码在MATLAB中，创建MIDAS（Mixed Data Sampling）模型需要进行几个步骤，包括定义MIDAS回归模型、生成MIDAS权重并执行模型拟合。

对于多元MIDAS和单变量MIDAS，虽然原理相同，但是在创建模型和权重生成时可能会有一些不同。

以下是一个单变量MIDAS模型和多元MIDAS模型的示例代码：注意：在使用这些代码之前，你需要先安装并导入Econometrics Toolbox，这个工具箱包含执行MIDAS分析的函数。

**单变量MIDAS模型**```matlab%导入数据data=readmatrix('your_data.csv');%替换为你的数据文件路径%定义MIDAS回归模型%假设我们有一个滞后因变量y和两个MIDAS滞后自变量x1和x2%在此例中，我们使用一个滞后阶数为4的MIDAS回归模型endog=y(:,1:4);exog=[ones(size(endog,1),1)x1(:,1:4)x2(:,1:4)];midas_model=midas_regress(endog,exog);%生成MIDAS权重%在此例中，我们使用一个滞后阶数为4的MIDAS权重生成器weights=midas_weights(endog,exog,4);%拟合模型results=estimate(midas_model,weights);```**多元MIDAS模型**对于多元MIDAS模型，你需要为每个因变量定义一个MIDAS回归模型，然后使用一个共同的MIDAS权重生成器。

假设你有一个包含两个因变量y1和y2的数据集：```matlab%导入数据data=readmatrix('your_data.csv');%替换为你的数据文件路径%定义MIDAS回归模型endog1=y1(:,1:4);exog1=[ones(size(endog1,1),1)x1(:,1:4)];%只使用x1作为自变量midas_model1=midas_regress(endog1,exog1); endog2=y2(:,1:4);exog2=[ones(size(endog2,1),1)x2(:,1:4)];%只使用x2作为自变量midas_model2=midas_regress(endog2,exog2);%生成MIDAS权重weights=midas_weights([endog1;endog2],[exog1; exog2],4);%使用一个共同的滞后阶数为4的MIDAS权重生成器%拟合模型results=estimate(midas_model1,weights);%使用第一个模型的拟合结果results2=estimate(midas_model2,weights);%使用第二个模型的拟合结果```以上代码示例只是一个简单的开始。

在Matlab中，我们可以使用以下的代码实现Grey-DEMATEL模型。

这段代码的主要目的是求解输入矩阵中的直接影响矩阵。

注意：这只是一个基础的示例，具体的代码可能需要根据实际的需求和问题进行修改。

```Matlab
function impactMatrix = grey_dematel(inputMatrix)
% 初始化参数
n = size(inputMatrix, 1);
impactMatrix = zeros(n);
% 计算直接影响矩阵
for i = 1:n
for j = 1:n
if inputMatrix(i, j) > 0
% 计算直接影响矩阵
impactMatrix(i, j) = inputMatrix(i, j) / sum(inputMatrix(i, :));
end
end
end
end
```
这个函数接受一个输入矩阵，然后计算并返回直接影响矩阵。

这个矩阵表示了每个因素对其他因素的影响程度。

注意：这个函数假设输入矩阵中的所有值都是非负的，因为在实际的Grey-DEMATEL模型中，我们只考虑非负的影响。

如果输入矩阵中有负值，那么这个函数可能无法正确地计算直接影响矩阵。

数学模型程序代码-M a t l a b-姜启源-第一章-建立数学模型-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN第1章 建立数学模型1.（求解，编程）如何施救药物中毒p10~11人体胃肠道和血液系统中的药量随时间变化的规律（模型）：d ,(0)1100d (,0)d ,(0)0d xx x ty x y y tλλμλμ⎧=-=⎪⎪>⎨⎪=-=⎪⎩ 其中，x (t )为t 时刻胃肠道中的药量，y (t )为t 时刻血液系统中的药量，t =0为服药时刻。

1.1（求解）模型求解p10~11要求：① 用MATLAB 求解微分方程函数dsolve 求解该微分方程（符号运算）。

② 用MATLAB 的化简函数simplify 化简所得结果。

③ 结果与教材P11上的内容比较。

提示：dsolve 和simplify 的用法可用help 查询。

建议在命令窗口中操作。

1.2（编程）结果分析p11已知λ=0.1386, μ=0.1155，将上题中得到x (t )和y (t )两条曲线画在同一个图形窗口内。

参考图形如下。

MATLAB命令plot, fplot, hold on/off, grid on/off, xlabel, ylabel, text 。

★ 编写的程序和运行结果：2.（编程，验证）商人们怎样安全过河p8~9三名商人各带一个随从乘船渡河，一只小船只能容纳二人，由他们自己划行。

随从们密约，在河的任一岸，一旦随从的人数比商人多，就杀人越货。

但是如何乘船的大权掌握在商人们手中。

商人们怎样才能安全渡河呢？[模型构成]决策：每一步（此岸到彼岸或彼岸到此岸）船上的人员。

要求：在安全的前提下（两岸的随从数不比商人多），经有限步使全体人员过河。

x k第k次渡河前此岸的商人数y k第k次渡河前此岸的随从数x k , y k=0,1,2,3; k=1,2,⋯过程的状态s k=(x k , y k)允许状态集合S={(x, y)|x=0, y=0,1,2,3; x=3, y=0,1,2,3; x=y=1,2}u k第k次渡船上的商人数v k第k次渡船上的随从数u k , v k=0,1,2; k=1,2,⋯决策d k=(u k , v k)允许决策集合D={(u , v)|u+v =1, 2}状态转移律s k+1=s k+(-1)k d k[多步决策问题]求d k∈D(k=1, 2, ⋯, n), 使s k∈S, 并按转移律由s1=(3,3) 到达s n+1=(0,0)。

function out=MMSmteam(s,m,mu1,mu2,T)%M/M/S/m排队模型%s——修理工个数%m——机器源数%T——时间终止点%mu1——机器离开-到达时间服从指数分布%mu2——修理时间服从指数分布%事件表：% p_s——修理工空闲概率% arrive_time——机器到达事件% leave_time——机器离开事件%mintime——事件表中的最近事件%current_time——当前时间%L——队长%tt——时间序列%LL——队长序列%c——机器到达时间序列%b——修理开始时间序列%e——机器离开时间序列%a_count——到达机器数%b_count——修理机器数%e_count——损失机器数%初始化arrive_time=exprnd(mu1,1,m);arrive_time=sort(arrive_time);leave_time=[];current_time=0;L=0;LL=[L];tt=[current_time];c=[];b=[];e=[];a_count=0;%循环while min([arrive_time,leave_time])<Tcurrent_time=min([arrive_time,leave_time]);tt=[tt,current_time]; %记录时间序列if current_time==min(arrive_time) %机器到达子过程arrive_time(1)=[]; % 从事件表中抹去机器到达事件a_count=a_count+1; %累加到达机器数if L<s %有空闲修理工L=L+1; %更新队长c=[c,current_time];%记录机器到达时间序列b=[b,current_time];%记录修理开始时间序列leave_time=[leave_time,current_time+exprnd(mu2)];%产生新的机器离开事件leave_time=sort(leave_time);%离开事件表排序else %无空闲修理工L=L+1; %更新队长c=[c,current_time];%记录机器到达时间序列endelse %机器离开子过程leave_time(1)=[];%从事件表中抹去机器离开事件arrive_time=[arrive_time,current_time+exprnd(mu1)];arrive_time=sort(arrive_time);%到达事件表排序e=[e,current_time];%记录机器离开时间序列if L>s %有机器等待L=L-1; %更新队长b=[b,current_time];%记录修理开始时间序列leave_time=[leave_time,current_time+exprnd(mu2)];%产生新的机器离开事件leave_time=sort(leave_time);%离开事件表排序else %无机器等待L=L-1; %更新队长endendLL=[LL,L]; %记录队长序列endWs=sum(e-c(1:length(e)))/length(e);Wq=sum(b-c(1:length(b)))/length(b);Wb=sum(e-b(1:length(e)))/length(e);Ls=sum(diff([tt,T]).*LL)/T;Lq=sum(diff([tt,T]).*max(LL-s,0))/T;p_s=1.0/(factorial(m)/factorial(m).*(mu2/mu1)^0+factorial(m)/factorial(m-1).*(mu2/mu1)^1+fact orial(m-2)/factorial(m-1).*(mu2/mu1)^2+factorial(m)/factorial(m-2).*(mu2/mu1)^2+factorial(m)/ factorial(m-4).*(mu2/mu1)^4+factorial(m)/factorial(m-5).*(mu2/mu1)^5);fprintf('修理工空闲概率:%d\n',p_s)%修理工空闲概率fprintf('到达机器数:%d\n',a_count)%到达机器数fprintf('平均逗留时间:%f\n',sum(e-c(1:length(e)))/length(e))%平均逗留时间fprintf('平均等待时间:%f\n',sum(b-c(1:length(b)))/length(b))%平均等待时间fprintf('平均修理时间:%f\n',sum(e-b(1:length(e)))/length(e))%平均修理时间fprintf('平均队长:%f\n',sum(diff([tt,T]).*LL)/T)%平均队长fprintf('平均等待队长:%f\n',sum(diff([tt,T]).*max(LL-s,0))/T)%平均等待队长for i=0:mp(i+1)=sum((LL==i).*diff([tt,T]))/T;%队长为i的概率fprintf('队长为%d的概率:%f\n',i,p(i+1));fprintf('机器不能马上得到修理的概率:%f\n',1-sum(p(1:s)))%机器不能马上得到修理的概率out=[Ws,Wq,Wb,Ls,Lq,p];function out=MMSkteam(s,k,mu1,mu2,T)%多服务台%s——服务台个数%k——最大顾客等待数%T——时间终止点%mu1——到达时间间隔服从指数分布%mu2——服务时间服从指数分布%事件表：% arrive_time——顾客到达事件% leave_time——顾客离开事件%mintime——事件表中的最近事件%current_time——当前时间%L——队长%tt——时间序列%LL——队长序列%c——顾客到达时间序列%b——服务开始时间序列%e——顾客离开时间序列%a_count——到达顾客数%b_count——服务顾客数%e_count——损失顾客数%初始化arrive_time=exprnd(mu1);leave_time=[];current_time=0;L=0;LL=[L];tt=[current_time];c=[];b=[];e=[];a_count=0;b_count=0;e_count=0;while min([arrive_time,leave_time])<Tcurrent_time=min([arrive_time,leave_time]);tt=[tt,current_time]; %记录时间序列if current_time==arrive_time %顾客到达子过程arrive_time=arrive_time+exprnd(mu1); % 刷新顾客到达事件a_count=a_count+1; %累加到达顾客数if L<s %有空闲服务台L=L+1; %更新队长b_count=b_count+1;%累加服务顾客数c=[c,current_time];%记录顾客到达时间序列b=[b,current_time];%记录服务开始时间序列leave_time=[leave_time,current_time+exprnd(mu2)];%产生新的顾客离开事件leave_time=sort(leave_time);%离开事件表排序elseif L<s+k %有空闲等待位L=L+1; %更新队长b_count=b_count+1;%累加服务顾客数c=[c,current_time];%记录顾客到达时间序列else %顾客损失e_count=e_count+1;%累加损失顾客数endelse %顾客离开子过程leave_time(1)=[];%从事件表中抹去顾客离开事件e=[e,current_time];%记录顾客离开时间序列if L>s %有顾客等待L=L-1; %更新队长b=[b,current_time];%记录服务开始时间序列leave_time=[leave_time,current_time+exprnd(mu2)];leave_time=sort(leave_time);%离开事件表排序else %无顾客等待L=L-1; %更新队长endendLL=[LL,L]; %记录队长序列endWs=sum(e-c(1:length(e)))/length(e);Wq=sum(b-c(1:length(b)))/length(b);Wb=sum(e-b(1:length(e)))/length(e);Ls=sum(diff([tt,T]).*LL)/T;Lq=sum(diff([tt,T]).*max(LL-s,0))/T;fprintf('到达顾客数:%d\n',a_count)%到达顾客数fprintf('服务顾客数:%d\n',b_count)%服务顾客数fprintf('损失顾客数:%d\n',e_count)%损失顾客数fprintf('平均逗留时间:%f\n',Ws)%平均逗留时间fprintf('平均等待时间:%f\n',Wq)%平均等待时间fprintf('平均服务时间:%f\n',Wb)%平均服务时间fprintf('平均队长:%f\n',Ls)%平均队长fprintf('平均等待队长:%f\n',Lq)%平均等待队长if k~=inffor i=0:s+kp(i+1)=sum((LL==i).*diff([tt,T]))/T;%队长为i的概率fprintf('队长为%d的概率:%f\n',i,p(i+1));endelsefor i=0:3*sp(i+1)=sum((LL==i).*diff([tt,T]))/T;%队长为i的概率fprintf('队长为%d的概率:%f\n',i,p(i+1));endendfprintf('顾客不能马上得到服务的概率:%f\n',1-sum(p(1:s)))%顾客不能马上得到服务的概率out=[Ws,Wq,Wb,Ls,Lq,p];。