【MATLAB第3期】源码分享#数学建模常用算法程序整理

- 本次分享内容包含神经网络、层次分析、移动平均、聚类、非线性优化问题、常微分方程问题、主成分分析、自动元胞机、图论、排队问题等。本次分享MATLAB及PYTHON编程语言解决数学建模问题使用的基本源代码，整理不易，求个关注。

**另赠送《MATLAB+神经网络43个案例分析》源代码&数据，关注加私聊“资源领取”，可回复下载链接：

《MATLAB 神经网络43个案例分析》目录
第1章 BP神经网络的数据分类——语音特征信号分类
第2章 BP神经网络的非线性系统建模——非线性函数拟合
第3章 遗传算法优化BP神经网络——非线性函数拟合
第4章 神经网络遗传算法函数极值寻优——非线性函数极值寻优
第5章 基于BP_Adaboost的强分类器设计——公司财务预警建模
第6章 PID神经元网络解耦控制算法——多变量系统控制
第7章 RBF网络的回归--非线性函数回归的实现
第8章 GRNN网络的预测----基于广义回归神经网络的货运量预测
第9章 离散Hopfield神经网络的联想记忆——数字识别
第10章 离散Hopfield神经网络的分类——高校科研能力评价
第11章 连续Hopfield神经网络的优化——旅行商问题优化计算
第12章 初始SVM分类与回归
第13章 LIBSVM参数实例详解
第14章 基于SVM的数据分类预测——意大利葡萄酒种类识别
第15章 SVM的参数优化——如何更好的提升分类器的性能
第16章 基于SVM的回归预测分析——上证指数开盘指数预测.
第17章 基于SVM的信息粒化时序回归预测——上证指数开盘指数变化趋势和变化空间预测
第18章 基于SVM的图像分割-真彩色图像分割
第19章 基于SVM的手写字体识别
第20章 LIBSVM-FarutoUltimate工具箱及GUI版本介绍与使用
第21章 自组织竞争网络在模式分类中的应用—患者癌症发病预测
第22章 SOM神经网络的数据分类--柴油机故障诊断
第23章 Elman神经网络的数据预测----电力负荷预测模型研究
第24章 概率神经网络的分类预测--基于PNN的变压器故障诊断
第25章 基于MIV的神经网络变量筛选----基于BP神经网络的变量筛选
第26章 LVQ神经网络的分类——乳腺肿瘤诊断
第27章 LVQ神经网络的预测——人脸朝向识别
第28章 决策树分类器的应用研究——乳腺癌诊断
第29章 极限学习机在回归拟合及分类问题中的应用研究——对比实验
第30章 基于随机森林思想的组合分类器设计——乳腺癌诊断
第31章 思维进化算法优化BP神经网络——非线性函数拟合
第32章 小波神经网络的时间序列预测——短时交通流量预测
第33章 模糊神经网络的预测算法——嘉陵江水质评价
第34章 广义神经网络的聚类算法——网络入侵聚类
第35章 粒子群优化算法的寻优算法——非线性函数极值寻优
第36章 遗传算法优化计算——建模自变量降维
第37章 基于灰色神经网络的预测算法研究——订单需求预测
第38章 基于Kohonen网络的聚类算法——网络入侵聚类
第39章 神经网络GUI的实现——基于GUI的神经网络拟合、模式识别、聚类
第40章 动态神经网络时间序列预测研究——基于MATLAB的NARX实现
第41章 定制神经网络的实现——神经网络的个性化建模与仿真
第42章 并行运算与神经网络——基于CPU/GPU的并行神经网络运算
第43章 神经网络高效编程技巧——基于MATLAB R2012b新版本特性的探讨

1.神经网络（MATLAB）【免费，源代码已验证】

%（1）普通BP神经网络

banana=[1 2 3 1 2 3 4 5 6 2 ];
pen=[2 3 7 2 3 4 5 3 9 2 ];
potato_actual=(banana*1+pen*2)/3;
potato_actual=[1.6 2.6 5.6 1.6 2.6 3.6 4.6 3.6 8.0 2.0];[input,ps1]=mapminmax([banana;pen]);
[target,ps2]=mapminmax([potato_actual]);net=newff(input,target,6,{'tansig','purelin'},'trainlm');
net.trainParam.epochs=1000;
net.trainParam.goal=0.00001;
LP.lr=0.000001;net=train(net,input,target);
banana1=[2 1 3 5 9 9];
pen1=[1 2 8 2 10 3];
potato_actual1=(banana1*1+pen1*2)/3;
input1=mapminmax('apply',[banana1;pen1],ps1);output1=net(input1);
prediction1=mapminmax('reverse',output1,ps2);set(0,'defaultfigurecolor','w')
figure
plot(potato_actual1,'*','color',[222 87 18]/255);hold on
plot(prediction1,'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction1'),title('预测其他数据')
xlabel('potato1'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)figure
output=net(input);
prediction=mapminmax('reverse',output,ps2);plot(potato_actual,'*','color',[29 131 8]/255);hold on
plot(prediction,'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('预测本身10个数据')
xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)

%（2）PSO-BP神经网络
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 该代码为基于PSO和BP网络的预测，保存为PSO.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 清空环境
clc
clear%读取数据
%load data input_train input_test output_train output_test%节点个数
inputnum=2;
hiddennum=3;
outputnum=1;%训练数据
input1=[1 2 3 1 2 3 4 5 6 2 ]; %输入(2个指标)
input2=[2 3 7 2 3 4 5 3 9 2 ];
output1=[1.6 2.6 5.6 1.6 2.6 3.6 4.6 3.6 8.0 2.0]; %输出
[inputn,inputps]=mapminmax([input1;input2]); %输入数据归一化
[outputn,outputps]=mapminmax(output1); %输出数据归一化
%测试数据
inputt1=[1 3 5];
inputt2=[2 7 3];
outputt1=[1.6 5.5 4.6];
inputtn=mapminmax('apply',[inputt1;inputt2],inputps); %输入数据归一化
%预测数据
in1=[1 3 5];
in2=[2 7 3];
inn=mapminmax('apply',[in1;in2],inputps); %输入数据归一化
%构建网络
net=newff(inputn,outputn,hiddennum);
% 参数初始化
%粒子群算法中的两个参数
c1 = 1.49445;
c2 = 1.49445;maxgen=10;   % 进化次数<----------------------------------------
sizepop=20;   %种群规模<---------------------------------------
wmax=0.9;
wmin=0.4;Vmax=1;
Vmin=-1;
popmax=5;
popmin=-5;
Dim=inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum;%% 产生初始粒子和速度
for i=1:sizepop%随机产生一个种群pop(i,:)=5*rands(1,Dim);    %初始种群vov(i,:)=rands(1,Dim);  %初始化速度%计算适应度fitness(i)=fun(pop(i,:),inputnum,hiddennum,outputnum,net,inputn,outputn,inputps,outputps);   %染色体的适应度
end% 个体极值和群体极值
[bestfitness bestindex]=min(fitness);
zbest=pop(bestindex,:);   %全局最佳
gbest=pop;    %个体最佳
fitnessgbest=fitness;   %个体最佳适应度值
fitnesszbest=bestfitness;   %全局最佳适应度值%% 迭代寻优
for i=1:maxgen%粒子位置和速度更新for j=1:sizepopw=wmax-(wmax-wmin)*j/maxgen;%速度更新%length(gbest(j,:));%length(pop(j,1:Dim))vov(j,:) = w*vov(j,:) + c1*rand*(gbest(j,:) - pop(j,1:Dim)) + c2*rand*(zbest - pop(j,1:Dim));vov(j,find(vov(j,:)>Vmax))=Vmax;vov(j,find(vov(j,:)<Vmin))=Vmin;%种群更新pop(j,1:Dim)=pop(j,1:Dim)+0.5*vov(j,:);pop(j,find(pop(j,1:Dim)>popmax))=popmax;pop(j,find(pop(j,1:Dim)<popmin))=popmin;%引入变异算子，重新初始化粒子if rand>0.9k=ceil(21*rand);pop(j,k)=rand;end%新粒子适应度值fitness(j)=fun(pop(j,1:Dim),inputnum,hiddennum,outputnum,net,inputn,outputn,inputps,outputps);end%%个体极值和群体极值更新for j=1:sizepop%个体最优更新if fitness(j) < fitnessgbest(j)gbest(j,:) = pop(j,1:Dim);fitnessgbest(j) = fitness(j);end%群体最优更新 if fitness(j) < fitnesszbestzbest = pop(j,1:Dim);fitnesszbest = fitness(j);endend%%每代最优值记录到yy数组中yy(i)=fitnesszbest;    end%% 结果分析
plot(yy)
title(['适应度曲线  ' '终止代数＝' num2str(maxgen)],'fontsize',12);
xlabel('进化代数','fontsize',12);ylabel('适应度','fontsize',12);x=zbest;
%% 把最优初始阀值权值赋予网络预测
% %用遗传算法优化的BP网络进行值预测
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=B2;%% BP网络训练
%网络进化参数
net.trainParam.epochs=100;
net.trainParam.lr=0.1;
net.trainParam.goal=0.0000001;%网络训练
[net,tr]=train(net,inputn,outputn);%% BP网络预测
%数据归一化
%预测训练数据
inputn_test=mapminmax('apply',[input1;input2],inputps);
an=sim(net,inputn_test);
%test_simu=mapminmax('reverse',[output1],outputps);
anss=mapminmax('reverse',[an],outputps);
error=output1-anss;figure(2)
plot(error)
title('仿真预测误差','fontsize',12);
xlabel('仿真次数','fontsize',12);ylabel('误差百分值','fontsize',12);plot(output1,'*','color',[29 131 8]/255);hold on
plot(anss,'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('预测本身数据')
xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)%预测测试数据
an=sim(net,inputtn);
anss=mapminmax('reverse',[an],outputps);
error=outputt1-anss;figure(3)
plot(error)
title('仿真预测误差','fontsize',12);
xlabel('仿真次数','fontsize',12);ylabel('误差百分值','fontsize',12);plot(outputt1,'*','color',[29 131 8]/255);hold on
plot(anss,'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('预测测试数据')
xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)%预测
an=sim(net,inn);
anss=mapminmax('reverse',[an],outputps)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%以下保存为fun.m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function error = fun(x,inputnum,hiddennum,outputnum,net,inputn,outputn,inputps,outputps)
%该函数用来计算适应度值
%x          input     个体
%inputnum   input     输入层节点数
%outputnum  input     隐含层节点数
%net        input     网络
%inputn     input     训练输入数据
%outputn    input     训练输出数据%error      output    个体适应度值%提取 BP神经网络初始权值和阈值，x为个体
w1=x(1:inputnum*hiddennum);
B1=x(inputnum*hiddennum+1:inputnum*hiddennum+hiddennum);
w2=x(inputnum*hiddennum+hiddennum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum);
B2=x(inputnum*hiddennum+hiddennum+hiddennum*outputnum+1:inputnum*hiddennum+hiddennum+hiddennum*outputnum+outputnum);%网络权值赋值
net.iw{1,1}=reshape(w1,hiddennum,inputnum);
net.lw{2,1}=reshape(w2,outputnum,hiddennum);
net.b{1}=reshape(B1,hiddennum,1);
net.b{2}=reshape(B2,outputnum,1);
%BP神经网络构建
net=newff(inputn,outputn,hiddennum);
net.trainParam.epochs=100;
net.trainParam.lr=0.1;
net.trainParam.goal=0.00001;
net.trainParam.show=100;
net.trainParam.showWindow=0;%BP神经网络训练
net=train(net,inputn,outputn);%网络训练
an=sim(net,inputn);
output=mapminmax('reverse',[outputn],outputps);
anss=mapminmax('reverse',[an],outputps);
error=sum(abs(anss-output));
end

2.层次分析法一致性检验（MATLAB）【免费，源代码已验证】

clc;
clear;
A = [1 5 3 71/5 1 1/3 31/3 3 1 51/7 1/3 1/5 1];
[m,n]=size(A);                     %获取指标个数
RI=[0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.51];
R=rank(A);                         %求判断矩阵的秩
[V,D]=eig(A);                      %求判断矩阵的特征值和特征向量，V特征值，D特征向量；
tz=max(D);
B=max(tz);                         %最大特征值
[row, col]=find(D==B);             %最大特征值所在位置
C=V(:,col);                        %对应特征向量
CI=(B-n)/(n-1);                    %计算一致性检验指标CI
CR=CI/RI(1,n);
if CR<0.10disp('CI=');disp(CI);disp('CR=');disp(CR);disp('对比矩阵A通过一致性检验，各向量权重向量Q为：');Q=zeros(n,1);for i=1:nQ(i,1)=C(i,1)/sum(C(:,1)); %特征向量标准化endQ                              %输出权重向量
elsedisp('对比矩阵A未通过一致性检验，需对对比矩阵A重新构造');
end
结果：
CI=107/2744
CR= 311/7178
对比矩阵A通过一致性检验，各向量权重向量Q为：
Q =     943/1669  484/4119  274/1045  91/1646  ```bash

3.移动平均（MATLAB）【免费，源代码已验证】

%（1）简单移动平均
% 简单移动平均只适合近期预测,并且使预测目标发展趋势变化不大的情况
clc,clear
y=[533.8 574.6 606.9 649.8 705.1 772.0 816.4 892.7 963.9 1015.1 1102.7];  % 已有的时间序列数据
m=length(y);
n=[2:5]; %n 为移动平均的项数
for i=1:length(n)
%由于 n 的取值不同，yhat 的长度不一致，下面使用了细胞数组 for j=1:m-n(i)+1yhat{i}(j)=sum(y(j:j+n(i)-1))/n(i); endpredict(i)=yhat{i}(end);error(i)=sqrt(mean((y(n(i)+1:m)-yhat{i}(1:end-1)).^2));
end
predict,errorfigure
plot(3:11,y(3:11),'-*','color',[29 131 8]/255);hold on
plot(3:11,yhat{1}(1:end-1),'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('移动平均项数为2时预测结果')
%xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)figure
plot(4:11,y(4:11),'-*','color',[29 131 8]/255);hold on
plot(4:11,yhat{2}(1:end-1),'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('移动平均项数为3时预测结果')
%xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)figure
plot(5:11,y(5:11),'-*','color',[29 131 8]/255);hold on
plot(5:11,yhat{3}(1:end-1),'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('移动平均项数为4时预测结果')
%xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)figure
plot(6:11,y(6:11),'-*','color',[29 131 8]/255);hold on
plot(6:11,yhat{4}(1:end-1),'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('移动平均项数为5时预测结果')
%xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)


%（2） 加权移动平均
% 加权移动平均，适合时间序列没有明显的趋势变动时的情况```bash
y=[6.35 6.20 6.22 6.66 7.15 7.89 8.72 8.94 9.28 9.8]; % 已有的时间序列数据
m=length(y);
w=[1/6;1/6;2/6;2/6]; % 权重
n=4; % 时间序列长度,必须有n==length(w)
yhat=[];
for i=1:m-n+1yhat(i)=y(i:i+n-1)*w;
end
err=abs(y(n+1:m)-yhat(1:end-1))./y(n+1:m)
T_err=1-sum(yhat(1:end-1))/sum(y(n+1:m))
predict=yhat(end)/(1-T_err) % 使用平均误差对结果进行修正
figure
plot(4:10,y(1,4:10),'-*','color',[29 131 8]/255);hold on
plot(4:10,yhat,'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','prediction')
title('时间序列为4时预测结果')
%xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)

%（3）二次移动平均
% 二次移动平均,适合时间序列出现直线增加或者减少的变动趋势
% 来自《数学建模算法与应用》P478
% 运用y的数据预测后两年

clc,clear
%load y.txt %把原始数据保存在纯文本文件 y.txt 中
y=[676 825 774 716 940 1159 1384 1524 1668 1688 1958 2031 2234 2566 2820 3006 3093 3277 3514 3770 4107]
m1=length(y);
n=6; %n 为移动平均的项数
for i=1:m1-n+1yhat1(i)=sum(y(i:i+n-1))/n;
end
yhat1
m2=length(yhat1);
for i=1:m2-n+1yhat2(i)=sum(yhat1(i:i+n-1))/n;
end
yhat2a21=2*yhat1(end)-yhat2(end)
b21=2*(yhat1(end)-yhat2(end))/(n-1)
predict1=a21+b21
predict2=a21+2*b21figure
plot(1:21,y,'-*','color',[29 131 8]/255);hold on
plot(6:21,yhat1,'-+','color','b');hold on
plot(11:21,yhat2,'-o','color',[244 208 0]/255,...
'linewidth',2,'MarkerSize',14,'MarkerEdgecolor',[138 151 123]/255);
legend('actua value','一次prediction','二次prediction')
title('移动项数为6时预测结果')
%xlabel('potato'),ylabel('weight')set(gca, 'Box', 'off', 'TickDir', 'out', 'TickLength', [.02 .02], ...'XMinorTick', 'on', 'YMinorTick', 'on', 'YGrid', 'on', ...'XColor', [.3 .3 .3], 'YColor', [.3 .3 .3],'LineWidth', 1)

4.聚类算法（MATLAB）【部分免费，源代码已验证】

%（1）简单的一维数据kmeans聚类

 A=[78    85    89    74    78    84    83];
B=[85    79    88    71    80    83    77];
C=[97    91    89    96    86    90    85];
D=[90    91    94    89    94    90    86];
E=[78    81    80    83    76    78    88];all=[A;B;C;D;E];
IDX=kmeans(all,2)  % 在这里分成两类
%在这里可以看到已经分成两类了，成绩较好的的IDX为1，稍逊的是2，

%结果：
%IDX =

  % 1       %1       %2       %2       % 1

%（2）二维数据kmeans聚类

x1=5*[randn(500,1)+5,randn(500,1)+5];
x2=5*[randn(500,1)+5,randn(500,1)-5];
x3=5*[randn(500,1)-5,randn(500,1)+5];
x4=5*[randn(500,1)-5,randn(500,1)-5];
x5=5*[randn(500,1),randn(500,1)];
all=[x1;x2;x3;x4;x5];  %生成2500*2的数组
plot(x1(:,1),x1(:,2),'r.');hold on
plot(x2(:,1),x2(:,2),'g.');...
plot(x3(:,1),x3(:,2),'k.');...
plot(x4(:,1),x4(:,2),'y.');...
plot(x5(:,1),x5(:,2),'b.');
IDX=kmeans(all,5); %这里用K-均值算法
y=pdist(all);
z=linkage(y);
t=cluster(z,'cutoff',1.2);
for k=1:2500text(all(k,1),all(k,2),num2str(IDX(k)));
end

%（3）R型聚类（对指标聚类）
%把下三角相关系数矩阵粘贴到纯文本文件ch.txt中 (即a数据)
%如
%1.0
%0.2 1.0
%0.4 0.5 1.0

clc,clear
a=[1.0 0 0;0.2 1.0 0;0.4 0.5 1.0];
d=1-abs(a);  %进行数据变换，把相关系数转化为距离
d=tril(d);   %提出d矩阵的下三角部分
b=nonzeros(d); %去掉d中的0元素
b=b';  %化成行向量
z=linkage(b,'complete');  %按最长距离法聚类
y=cluster(z,'maxclust',2);%把变量划分成两类，注:也可3类，底下记得修改
ind1=find(y==1);  %显示第一类对应的变量编号
ind2=find(y==2);  %显示第二类对应的变量编号
ind1=ind1';
ind2=ind2';
h=dendrogram(z);  %画聚类图
ind1,ind2
set(h,'Color','k','LineWidth',2.0);%把聚类图线的颜色修改成黑色，线宽加粗

%（4）Q型聚类（对样本聚类）

clc,clear;
a=[1,0;1,1;3,2;4,3;2,5];  %共5个样本,每个样本两个指标。欲将5个样本分为3类
a=zscore(a);              %数据标准化处理
y=pdist(a,'cityblock');   %求a的行向量之间的绝对距离
% http://blog.sciencenet.cn/blog-531885-589056.html 各种距离
yc=squareform(y);         %变换成距离方阵
z=linkage(y);             %产生等级聚类图
[h,t]=dendrogram(z);      %画聚类图
T=cluster(z,'maxclust',3);%把对象分成3份，参数可自行修改成2,4,5等，记得将下一行i值修改
for i=1:3;  tm=find(T==i);        %求第i类的对象  tm=reshape(tm,1,length(tm));%变成行向量  fprintf('第%d类的有%s\n',i,int2str(tm));%显示分类结果
end

第1类的有1 2第2类的有3 4第3类的有5
%（5）二维聚类（模糊聚类）

%% 1、随机产生初始聚类中心
% clc;
clear
%% 加载数据
load X
figure
plot(X(:,1),X(:,2),'o')
xlabel('横坐标X');ylabel('纵坐标Y');title('样本数据')
hold on
%% 进行模糊C均值聚类
% 设置幂指数为3，最大迭代次数为20，目标函数的终止容限为1e-6
options=[3,20,1e-6,0];
% 调用fcm函数进行模糊C均值聚类，返回类中心坐标矩阵center，隶属度矩阵U，目标函数值obj_fcn
cn=4; %聚类数
[center,U,obj_fcn]=fcm(X,cn,options);
Jb=obj_fcn(end)
maxU = max(U);
index1 = find(U(1,:) == maxU)
index2 = find(U(2, :) == maxU)
index3 = find(U(3, :) == maxU)
%% 分类情况
% 在前三类样本数据中分别画上不同记号 不加记号的就是第四类了
line(X(index1,1), X(index1, 2), 'linestyle', 'none', 'marker', 'x', 'color', 'g');
line(X(index2,1), X(index2, 2), 'linestyle', 'none', 'marker', '*', 'color', 'r');
line(X(index3,1), X(index3, 2), 'linestyle', 'none', 'marker', '+', 'color', 'b');
%% 画出聚类中心
plot(center(:,1),center(:,2),'v')
hold off

```![在这里插入图片描述](https://img-blog.csdnimg.cn/729093c79e154aee8007b20194666ac3.png)%（6）二维聚类（模糊&遗传&退火聚类）【源代码已验证】```bash
clc
load X
m=size(X,2);% 样本特征维数
% 中心点范围[lb;ub]
lb=min(X);
ub=max(X);
%% 模糊C均值聚类参数
% 设置幂指数为3，最大迭代次数为20，目标函数的终止容限为1e-6
options=[3,20,1e-6];
% 类别数cn
cn=4;
%% 模拟退火算法参数
q =0.8;     % 冷却系数
T0=100;    % 初始温度
Tend=99.999;  % 终止温度
%% 定义遗传算法参数
sizepop=10;               %个体数目(Numbe of individuals)
MAXGEN=100;                %最大遗传代数(Maximum number of generations)
NVAR=m*cn;                %变量的维数
PRECI=10;                 %变量的二进制位数(Precision of variables)
pc=0.7;
pm=0.01;

5.非线性优化问题（MATLAB）【免费，源代码已验证】

%（1）非线性优化问题（非线性表达式 s.t. 约束条件）

% 非线性优化问题（非线性表达式 s.t. 约束条件）
% 例如求非线性表达式为 f(x)=\sum{i=1}^{11} b_i * \log(x_i+t_i) - b_i * \log(t_i)的最大值
% 约束条件为x_i>=0 (i=1,2,...11) 以及 \sum_{i=1}^{11} x_i=100000000% gaoptimset('PlotFcns',@gaplotbestf)为options的参数，可以画出下降图，如果画出来不合适，可以使用
% gaoptimset('Generations',100,'PlotFcns',@gaplotbestf)来限定迭代次数，这样输出的图的横坐标就限制在100内
% 在工具箱里面也可以设置plot的参数% 如果有非线性约束，见http://blog.csdn.net/qq_37043191/article/details/77898335
% 如果有整数约束，则在GUI的IntConIndex中填写如[1 3 5]（表示x1,x3,x5为整数）
% 01问题则在IntConIndex中填写[1 2 ....]，再将每个变量上限下限填为1和0% 以下内容保存为fun.m
function f =fun(x)b= [0.903, 0.501, 0.655 , 1.075 , 0.965 , 0.629 , 1.429 , 1.275 , 0.755 , 1.701 , 1.219];t= [60202000,22027551,105132543,2929453717 ,71205991,95780211,36908000,48132000,111949554,62312533,51877376];f=b(1)*log(t(1)+x(1))-b(1)*log(t(1));for i=2:11f = f+b(i)*log(t(i)+x(i))-b(i)*log(t(i));  endf=-f; % 转化为求最小值
end% [x,fval] =fmincon('函数名',ga(@函数名, 变量个数, A, b, Aeq, beq),A ,b ,Aeq ,beq);
% [x,fval] =ga(@fun,11,-eye(11,11),zeros(11,1),ones(1,11),[100000000],[],[],[],gaoptimset('PlotFcns',@gaplotbestf));
% [x,fval] =fmincon('fun',ga(@fun,11,-eye(11,11),zeros(11,1),ones(1,11),[100000000]),-eye(11,11),zeros(11,1),ones(1,11),[100000000]);

%（2）非线性方程无约束优化(遗传算法&fminunc)

% 待优化方程保存为fun.m
function f=fun(x)f=((33*x(3))/10 - (33*x(2))/10 + 231/5)*(x(3) - x(2)+ 14) + ((37*x(2))/10 - (37*x(1))/10 + 851/10)*(x(2) - x(1) + 23) + (29*x(1)^2)/10 +(5*x(2)^2)/2 + 4*x(3)^2;
end% [x,fval] =fminunc('fun',ga(@fun,3))

%（3）非线性方程无约束优化(模拟退火&fminunc)

function sol% 保存为fun.mfunction f=fun(x)f=((33*x(3))/10 - (33*x(2))/10 + 231/5)*(x(3) - x(2)+ 14) + ((37*x(2))/10 - (37*x(1))/10 + 851/10)*(x(2) - x(1) + 23) + (29*x(1)^2)/10 +(5*x(2)^2)/2 + 4*x(3)^2;end% [x,fval]=simulannealbnd(@fun,[0;0;0])[x,fval] =fminunc('fun',simulannealbnd(@fun,[0;0;0]));
end

%（4）非线性方程非线性约束优化

% 求解非线性约束问题
% 待优化表达式f(x, y)=x^3-y^3+x*y+2*x^2
% 约束条件为x^2+y^2<=6和x*y=2
% 创建以下函数并保存为fun.m
function f=fun(x)f=x(1)^3-x(2)^3+x(1)*x(2)+2*x(1)^2;
end% 创建约束条件的函数文件并保存为fcontr.m
function [c, d]=fcontr(x)c=x(1)^2+x(2)^2-6; % c为不等式约束，需转化为<=0的格式d=x*y-2; % d为等式约束，需化为=0的格式
end% [x fval exitflag]=fmincon('fun',[1 2],[],[],[],[],[],[],'fcontr');

%（5）非线性方程无约束（进化算法） PYTHON
含WOA / PSO / MFO / MVO / GWO /FFA /CS /BAT等进化算法。
整理中，后续会在第三方平台免费分享。
1.PSO论文参考《Particle Swarm Optimization》
百度学术：
2.GWO论文参考《Grey Wolf Optimizer》：“the GWO algorithm is able to provide very competitive results compared to well-known meta-heuristics.”
（1）来源：
（2）google学术：
3.EvoloPy论文参考《An Open-Source Nature-Inspired Optimization Framework in Python》：“EvoloPy is an open source and cross-platform Python framework that implements a wide range of classical and recent nature-inspired metaheuristic algorithms.”
（1）百度学术：

6.解微分方程（MATLAB）【免费，源代码已验证】

%（1）一阶常微分方程组(ode45)

function sol% dx/dt=y% dy/dt=z% dz/dt=x^3-x-y-0.3175z% 初值x(0)=0,y(0)=1,z(0)=2function dx=fun(t,x)dx(1)=x(2);dx(2)=x(3);dx(3)=x(1)^3-x(1)-x(2)-0.3175*x(3);dx=dx(:);end% [t,x]=ode45(@fun,求解区间,初值)[t,x]=ode45(@fun,[0:0.01:2],[0,1,2])% plot3(x(:,1),x(:,2),x(:,3))
end

%（2）一阶常微分方程(ode45)

function sol% dy=(y+3t)/t^2% 初值y(0)=2function dx=fun(t,x)dx(1)=(x(1)+3*t)/t^2;dx=dx(:);end[t,x]=ode45(@fun,[1:0.01:4],[-2])plot(t,x)
end

%（3）二阶常微分方程(ode45)

function sol% 高阶常微分方程 d2y=-t*y+exp(t)*dy+3*sin(2*t)% 初值 y(0)=2,dy(0)=8tspan=[3.9 4.0]; %求解区间y0=[2 8]; %初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y'' ''=-t*y + e^t*y'' +3sin2t')xlabel('t')ylabel('y')function y=odefun(t,x)y=zeros(2,1); % 列向量 x(1)代表y,x(2)代表dyy(1)=x(2);    % 这里y(1)代表dy,y(2)代表d2yy(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);end
end

7.主成分分析（PYTHON）【免费，源代码已验证】

#coding=utf-8
from numpy import *'''通过方差的百分比来计算将数据降到多少维是比较合适的，
函数传入的参数是特征值和百分比percentage，返回需要降到的维度数num'''
def eigValPct(eigVals,percentage):sortArray=sort(eigVals) #使用numpy中的sort()对特征值按照从小到大排序sortArray=sortArray[-1::-1] #特征值从大到小排序arraySum=sum(sortArray) #数据全部的方差arraySumtempSum=0num=0for i in sortArray:tempSum+=inum+=1if tempSum>=arraySum*percentage:return num

'''pca函数有两个参数，其中dataMat是已经转换成矩阵matrix形式的数据集，列表示特征；
其中的percentage表示取前多少个特征需要达到的方差占比，默认为0.9'''
def pca(dataMat,percentage=0.9):meanVals=mean(dataMat,axis=0)  #对每一列求平均值，因为协方差的计算中需要减去均值meanRemoved=dataMat-meanValscovMat=cov(meanRemoved,rowvar=0)  #cov()计算方差eigVals,eigVects=linalg.eig(mat(covMat))  #利用numpy中寻找特征值和特征向量的模块linalg中的eig()方法k=eigValPct(eigVals,percentage) #要达到方差的百分比percentage，需要前k个向量eigValInd=argsort(eigVals)  #对特征值eigVals从小到大排序eigValInd=eigValInd[:-(k+1):-1] #从排好序的特征值，从后往前取k个，这样就实现了特征值的从大到小排列redEigVects=eigVects[:,eigValInd]   #返回排序后特征值对应的特征向量redEigVects（主成分）lowDDataMat=meanRemoved*redEigVects #将原始数据投影到主成分上得到新的低维数据lowDDataMatreconMat=(lowDDataMat*redEigVects.T)+meanVals   #得到重构数据reconMatreturn lowDDataMat,reconMatdata=array([[2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1],[2.4,0.7,2.9,2.2,3,2.7,1.6,1.1,1.6,0.9]])
l, r= pca(data.T)
print("降维后的结果：")
print(l)
print("用降维后的结果重构的：")
print(r)
8.自动元胞机（MATLAB） 【源代码已验证】
%（1）多车道元胞机（可换道、随机变慢）
data=[];
%%%%%%%%%%%%%%%%%%%set%%%%%%%%%%%%%%%%
length=100;
lane_num=8;
p_new_car=0.8;
p_auto_car=0;
max_speed=4;ca_show=zeros(length,100,lane_num); %第一维为车道长度,第二维为时间维,第三维为车道
for j=1:100[loc_matrix,v,vmax]=new_cars(zeros(length,lane_num),zeros(length,lane_num),zeros(length,lane_num),p_new_car,p_auto_car,max_speed); % 注意这里改了在move_forward里面也要改traffic_count=0;for i=1:1000[loc_matrix,v,vmax,traffic_count]=move_forward(loc_matrix,v,vmax,100,traffic_count);ca_show(:,2:end,:)=ca_show(:,1:end-1,:);for i=1:lane_numca_show(:,1,i)=loc_matrix(:,i);end%imshow([loc_matrix(:,2) 0.5*ones(length(loc_matrix(:,1)),1) loc_matrix(:,3) 0.5*ones(length(loc_matrix(:,1)),1) loc_matrix(:,4) 0.5*ones(length(loc_matrix(:,1)),1)...%    loc_matrix(:,5) 0.5*ones(length(loc_matrix(:,1)),1) loc_matrix(:,6) 0.5*ones(length(loc_matrix(:,1)),1) loc_matrix(:,7)],'InitialMagnification','fit')%figure(1)imshow(ca_show(:,:,1),'InitialMagnification','fit')xlabel('t')ylabel('one lane')%figure(2)%imshow(loc_matrix,'InitialMagnification','fit')pause(0.03)enddata=[data traffic_count];
end
disp(mean(data))

%（2）NS规则下模拟车流量随车密度的变化（单环）

clc;
clear all;
% 求车辆的流量随车辆密度(从0到1)变化的规律
T=3030;  %3030
P=0.3; % 随机慢化的概率
v_max=6; % vn在0,1,...6中取值
L=2000; % 单车道,道路总长为2000个格子(是个环,头尾相接),放N个车
dens=0.002; %给定初始车辆密度
p=1; %统计流量密度数组
while dens<=1 N=fix(dens*L); %车辆数目,N一定小于Lm=1;% 产生N辆车的初始随机速度v_matrix=randperm(N); % randperm产生1,2,...,N的打乱了的数字for i=1:Nv_matrix(i)=mod(v_matrix(i),v_max+1); % 将速度投射到合理范围end% 产生初始随机位置,N个车放入单车道中(L个格子)[a,b]=find(randperm(L)<=N);loc_matrix=b; % 位置从小到大排列,如5辆车时,434,509,1513,1917%变化规则 for i=1:T%定义车头间距,如果两辆车在相邻格子中,则车头间距认为是0if loc_matrix(N)>loc_matrix(1) % (以1到2000为序),计算1917格子的车-->2000-->1-->434格子车头的距离% 这里的N的一开始的编号,之后N对应的位置并不一定是1-2000中最大的headways(N)=L-loc_matrix(N)+loc_matrix(1)-1;elseheadways(N)=loc_matrix(1)-loc_matrix(N)-1;endfor j=N-1:-1:1if loc_matrix(j+1)>loc_matrix(j)headways(j)=loc_matrix(j+1)-loc_matrix(j)-1;elseheadways(j)=L+loc_matrix(j+1)-loc_matrix(j)-1;endend%速度变化 v_matrixNS1=min([v_max-1,v_matrix(1),max(0,headways(1)-1)]);%随机慢化概率 %最大速度 %网格的数量v_matrix(N)=min([v_max,v_matrix(N)+1,headways(N)+v_matrixNS1]); %NS规则下第N辆车的速度估计值for j=N-1:-1:1 v_matrixNS=min([v_max-1,v_matrix(j+1),max(0,headways(j+1)-1)]); %NS规则下前一辆车的速度估计值v_matrix(j)=min([v_max,v_matrix(j)+1,headways(j)+v_matrixNS]); %NS规则下第j辆车的前一辆车的速度估计值end%以概率P随机慢化if rand()<P;v_matrix=max(v_matrix-1,0); %随机慢化规则下的速度变化end%位置更新for j=N:-1:1loc_matrix(j)=loc_matrix(j)+v_matrix(j);if loc_matrix(j)>=Lloc_matrix(j)=loc_matrix(j)-L; %NS规则下第j辆车的位置更新endend%采集数据作图if i>L+1000 %采用每组的后30个变量取平均speed(m)=sum(v_matrix)/N; %求取平均速度m=m+1;endendflow(p)=(sum(speed)/30)*dens; %不同密度下的流量数组 density(p)=dens;dens=dens+0.01;p=p+1;
end
plot(density,flow)

%（3）多车道从左超车模型（MATLAB）【源代码已验证】

clc;
clear all;
close all;B=4;             %The number of the lanes
plazalength=50;  %The length of the simulating highways
h=NaN;           %h is the handle of the image[plaza,v]=create_plaza(B,plazalength);
h=show_plaza(plaza,h,0.1);iterations=1000;    % 迭代次数
probc=0.1;          % 车辆的密度
probv=[0.1 1];      % 两种车流的密度分布
probslow=0.3;       % 随机慢化的概率
Dsafe=1;            % 表示换道事车至少与后面车距离多少个单位才算安全
VTypes=[1,2];       %道路上一共有几种最大速度不同的车辆,速度是什么
[plaza,v,vmax]=new_cars(plaza,v,probc,probv,VTypes);%一开始就在车道上布置车辆，做周期循环驾驶，也方便观察流量密度之间的关系
```![在这里插入图片描述](https://img-blog.csdnimg.cn/66cd4925345a4045a2a519766432b368.png)%（4）自动驾驶汽车（MATLAB） 【源代码已验证】```bash
function Zmain()
%nl:车道长度;nc:车道数目;fp:车道入口处新进入车辆的概率;dt:仿真步长时间;nt:仿真时间;tcinfo:车辆行驶相关参数
[nl,nc,fp,dt,nt] = ZGetBasicInfo();
tcinfo=ZGetTrafficiInfo(); %随机慢化概率、元胞最大速度、随机换道概率、收费站位置、自动驾驶汽车比例
%生成元胞空间
cellspace=ZGenerateCellSpace(nl,nc,tcinfo);
%开始仿真
ZTrafficSimulating(cellspace,fp,dt,nt,tcinfo);
end

9.图论算法（MATLAB）【源代码已验证】

%（1）基于蚁群算法的三维路径规划算法（MATLAB）【源代码已验证】

%% 数据初始化%下载数据
load  HeightData HeightData%网格划分
LevelGrid=10;
PortGrid=21;%起点终点网格点
starty=10;starth=4;
endy=8;endh=5;
m=1;
%算法参数
PopNumber=10;         %种群个数
BestFitness=[];    %最佳个体%初始信息素
pheromone=ones(21,21,21);

%（2）蚁群算法的优化计算——旅行商问题(TSP)优化（MATLAB）【源代码已验证】

%% 清空环境变量
clear all
clc%% 导入数据
load citys_data.mat%% 计算城市间相互距离
n = size(citys,1);
D = zeros(n,n);
for i = 1:nfor j = 1:nif i ~= jD(i,j) = sqrt(sum((citys(i,:) - citys(j,:)).^2));elseD(i,j) = 1e-4;      endend
end%% 初始化参数
m = 50;                              % 蚂蚁数量
alpha = 1;                           % 信息素重要程度因子
beta = 5;                            % 启发函数重要程度因子
rho = 0.1;                           % 信息素挥发因子
Q = 1;                               % 常系数
Eta = 1./D;                          % 启发函数
Tau = ones(n,n);                     % 信息素矩阵
Table = zeros(m,n);                  % 路径记录表
iter = 1;                            % 迭代次数初值
iter_max = 200;                      % 最大迭代次数
Route_best = zeros(iter_max,n);      % 各代最佳路径
Length_best = zeros(iter_max,1);     % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);      % 各代路径的平均长度

最短距离:15601.9195
最短路径:14 12 13 11 23 16 5 6 7 2 4 8 9 10 3 18 17 19 24 25 20 21 22 26 28 27 30 31 29 1 15 14

%（3）免疫优化算法解决物流配送选址（给定需求坐标以及每个坐标需求量）（MATLAB）【源代码已验证】

%% 清空环境
clc
clear%% 算法基本参数
sizepop=50;           % 种群规模
overbest=10;          % 记忆库容量
MAXGEN=100;            % 迭代次数
pcross=0.5;           % 交叉概率
pmutation=0.4;        % 变异概率
ps=0.95;              % 多样性评价参数
length=6;             % 配送中心数
M=sizepop+overbest;%% step1 识别抗原,将种群信息定义为一个结构体
individuals = struct('fitness',zeros(1,M), 'concentration',zeros(1,M),'excellence',zeros(1,M),'chrom',[]);
%% step2 产生初始抗体群
individuals.chrom = popinit(M,length);
trace=[]; %记录每代最个体优适应度和平均适应度

%（4）PSO解TSP（给定城市坐标）（MATLAB）【源代码已验证】

%% 下载数据
data=load('eil51.txt');
cityCoor=[data(:,2) data(:,3)];%城市坐标矩阵,第一维是编号figure
plot(cityCoor(:,1),cityCoor(:,2),'ms','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g')
legend('城市位置')
ylim([4 78])
title('城市分布图','fontsize',12)
xlabel('km','fontsize',12)
ylabel('km','fontsize',12)

10.排队问题（MATLAB / PYTHON / C）【源代码已验证】

PYTHON及C相关的代码会在第三方平台免费分享下载。有需求者也可以关注我，私信。
%（1）多通道等待模拟（MATLAB ）【源代码已验证】

 [num,pass]=computing(1000)
csvwrite('data.csv',pass)
function [num,pass]=computing(tim0)seat=[0 0 0];%服务员属性pass=rand(1,4);%顾客信息：序号、到达时间、特殊要求时间、正常理发时间pass(5)=0;%理发员pass(6)=0;%离开时间pass(7)=0;%等待时间num=1;%服务人数初始化tim=0;%时间计数器temp=0;%顾客到达时间间隔

%（2）单通道等待模拟（MATLAB ）【免费，源代码已验证】

% 单通道等待模拟
% 货车夜间到达，白天卸货，每天只能卸货2车，若一天内到达数超过2车，那么就推迟到次日卸货
% 到达车数 0    1    2   3    4    5   ≥6
% 概 率  0.23 0.30 0.30 0.1 0.05 0.02 0.00
% 到达车数不服从泊松分布，服务时间也不服 从指数分布(这是定长服务时间)
% a1 表示产生的随机数，a2 表示到达的车数，a3 表示需要卸货车数，a4 表示实际卸货车数，a5 表示推迟卸货车数
clear
rand('state',sum(100*clock));
n=50000;
m=2;
a1=rand(n,1);
a2=a1; %a2初始化
a2(find(a1<0.23))=0;
a2(find(0.23<=a1&a1<0.53))=1;
a2(find(0.53<=a1&a1<0.83))=2;
a2(find(0.83<=a1&a1<0.93),1)=3;
a2(find(0.93<=a1&a1<0.98),1)=4;
a2(find(a1>=0.98))=5;
a3=zeros(n,1);
a4=zeros(n,1);
a5=zeros(n,1); %a2初始化
a3(1)=a2(1);
if a3(1)<=ma4(1)=a3(1);a5(1)=0;
elsea4(1)=m;a5(1)=a2(1)-m;
end
for i=2:na3(i)=a2(i)+a5(i-1);if a3(i)<=ma4(i)=a3(i);a5(i)=0;elsea4(i)=m;a5(i)=a3(i)-m;end
end
a=[a1,a2,a3,a4,a5];
sum(a)/n

%（3）matlab 排队含GUI界面（MATLAB ）【源代码已验证】

%   g为顾客到达时间服从的泊松分布的参数
%   h为顾客服务时间服从的负指数分布的参数
%   N为服务台数
%   Tmax为仿真时长
if get(handles.start,'String')=='退出'close;
elseset(handles.start,'String','退出');global N;global g;global h;global Tmax;%vf存放仿真速度global vf;%A为顾客的详细信息，其中按行依次为到达时间、服务时长、排队时长、服务台号,离开时间和进入服务时间global A;A=zeros(6,1);%离开时间排序global leave;leave=0;%S为服务台的信息，其中行代表服务台号，列分别代表状态（0为空闲，1为工作），和状态改变的时间global S;S=zeros(N,1);%s表示对应的服务台的状态改变数global s;s=zeros(1,N);%Se表示总的服务利用率，Sm表示表示服务台处于工作状态的台数global Se;Se=0;global Sm;Sm=0;%L为队列信息，行分别表示状态改变时间和相应的队长,Te表示平均等待时间，Le表示平均队长global L;L=[0;0];global Te;Te=0;global Le;Le=0;%t存放状态改变的时刻，顾客号，及相应的状态（1到达，0离开,2进入服务）global t;

0s时，1号顾客到达商场
0s时，1号顾客到1号服务台接受服务
4s时，2号顾客到达商场
4s时，3号顾客到达商场
4.3916s时，1号顾客从1号服务台离开
4.3916s时，2号顾客到1号服务台接受服务
4.879s时，2号顾客从1号服务台离开
4.879s时，3号顾客到1号服务台接受服务
5.0236s时，3号顾客从1号服务台离开
7s时，4号顾客到达商场