1.层次分析；
2.多属性决策；
3.图论Dijstra;
4.图论Floryd;
5.退火算法；
6.种群竞争；
1.层次分析法：
实现代码：
disp('请输入判断矩阵A(n阶)');
A=input('A=');
[n,n]=size(A);
x=ones(n,100);
y=ones(n,100);
m=zeros(1,100);
m(1)=max(x(:,1));
y(:,1)=x(:,1);
x(:,2)=A*y(:,1);
m(2)=max(x(:,2));
y(:,2)=x(:,2)/m(2);
p=0.0001;i=2;k=abs(m(2)-m(1));
while k>p
i=i+1;
x(:,i)=A*y(:,i-1);
m(i)=max(x(:,i));
y(:,i)=x(:,i)/m(i);
k=abs(m(i)-m(i-1));
end
a=sum(y(:,i));
w=y(:,i)/a;
t=m(i);
disp(w);
%以下是一致性检验
CI=(t-n)/(n-1);RI=[0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];
CR=CI/RI(n);
if CR<0.10
disp('此矩阵的一致性可以接受!');
disp('CI=');disp(CI);
disp('CR=');disp(CR);
end
测试数据：
[1, 1/2, 4, 3, 3;
2, 1, 7, 5, 5;
1/4, 1/7, 1, 1/2, 1/3;
1/3, 1/5, 2, 1, 1;
1/3, 1/5, 3, 1, 1;]
[1,2,5;
1/2,1,2;
1/5,1/2,1;]
[1,1/3,1/8;
3,1,1/3;
8,3,1;]
[1,1,3;
1,1,3;
1/3,1/3,1;]
[1,3,4;
1/3,1,1;
1/4,1,1;]
[1,1,1/4;
1,1,1/4;
4,4,1;]
2.多属性决策：
代码实现：
disp('请输入判断矩阵A(n阶)');
A=input('A=');
[n,n]=size(A);
x=ones(n,100);
y=ones(n,100);
m=zeros(1,100);
m(1)=max(x(:,1));
y(:,1)=x(:,1);
x(:,2)=A*y(:,1);
m(2)=max(x(:,2));
y(:,2)=x(:,2)/m(2);
p=0.0001;i=2;k=abs(m(2)-m(1));
while k>p
i=i+1;
x(:,i)=A*y(:,i-1);
m(i)=max(x(:,i));
y(:,i)=x(:,i)/m(i);
k=abs(m(i)-m(i-1));
end
a=sum(y(:,i));
w=y(:,i)/a;
t=m(i);
disp(w);
%以下是一致性检验
CI=(t-n)/(n-1);RI=[0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];
CR=CI/RI(n);
if CR<0.10
disp('此矩阵的一致性可以接受!');
disp('CI=');disp(CI);
disp('CR=');disp(CR);
end
测试数据：
[1 3 3 3 3;
1/3 1 1 1 1;
1/3 1 1 1 1;
1/3 1 1 1 1;
1/3 1 1 1 1;]3.灰色预测：
function []=greymodel(y)
% 本程序主要用来计算根据灰色理论建立的模型的预测值。
% 应用的数学模型是 GM(1,1)。
% 原始数据的处理方法是一次累加法。
y=input('请输入数据 ');
n=length(y);
yy=ones(n,1);
yy(1)=y(1);
for i=2:n
yy(i)=yy(i-1)+y(i);
end
B=ones(n-1,2);
for i=1:(n-1)
B(i,1)=-(yy(i)+yy(i+1))/2;
B(i,2)=1;
end
BT=B';
for j=1:n-1
YN(j)=y(j+1);
end
YN=YN';
A=inv(BT*B)*BT*YN;
a=A(1);
u=A(2);
t=u/a;
i=1:n+2;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;
yys(1)=y(1);
for j=n+2:-1:2
ys(j)=yys(j)-yys(j-1);
end
x=1:n;
xs=2:n+2;
yn=ys(2:n+2);
plot(x,y,'^r',xs,yn,'*-b');
det=0;
sum1=0;
sumpe=0;
for i=1:n
sumpe=sumpe+y(i);
end
pe=sumpe/n;
for i=1:n;
sum1=sum1+(y(i)-pe).^2;
end
s1=sqrt(sum1/n);
sumce=0;
for i=2:n
sumce=sumce+(y(i)-yn(i));
end
ce=sumce/(n-1);
sum2=0;
for i=2:n;
sum2=sum2+(y(i)-yn(i)-ce).^2;
end
s2=sqrt(sum2/(n-1));
c=(s2)/(s1);
disp(['后验差比值为：',num2str(c)]);
if c<0.35
disp('系统预测精度好')
else if c<0.5
disp('系统预测精度合格')
else if c<0.65
disp('系统预测精度勉强')
else
disp('系统预测精度不合格')
end
end
end
disp(['下个拟合值为 ',num2str(ys(n+1))]);
disp(['再下个拟合值为',num2str(ys(n+2))]);测试数据：
[724.57, 746.62, 778.27, 800.8, 827.75,871.1, 912.37, 954.28, 995.01, 1037.2]
[2.874,3.278,3.337,3.390,3.679]4.图论_dijstra算法：
tulun1.m
weight= [0 2 8 1 Inf Inf Inf Inf Inf Inf Inf;
2 0 6 Inf 1 Inf Inf Inf Inf Inf Inf;
8 6 0 7 5 1 2 Inf Inf Inf Inf;
1 Inf 7 0 Inf Inf 9 Inf Inf Inf Inf;
Inf 1 5 Inf 0 3 Inf 2 9 Inf Inf;
Inf Inf 1 Inf 3 0 4 Inf 6 Inf Inf;
Inf Inf 2 9 Inf 4 0 Inf 3 1 Inf;
Inf Inf Inf Inf 2 Inf Inf 0 7 Inf 9;
Inf Inf Inf Inf 9 6 3 7 0 1 2;
Inf Inf Inf Inf Inf Inf 1 Inf 1 0 4;
Inf Inf Inf Inf Inf Inf Inf 9 2 4 0;];
[dis, path]=dijkstra(weight,1, 11)
dijkstra.m
function [min,path]=dijkstra(w,start,terminal)
n=size(w,1); label(start)=0; f(start)=start;
for i=1:n
if i~=start
label(i)=inf;
end, end
s(1)=start; u=start;
while length(s)<n
for i=1:n
ins=0;
for j=1:length(s)
if i==s(j)
ins=1;
end,
end
if ins==0
v=i;
if label(v)>(label(u)+w(u,v))
label(v)=(label(u)+w(u,v));
f(v)=u;
end,
end,
end
v1=0;
k=inf;
for i=1:n
ins=0;
for j=1:length(s)
if i==s(j)
ins=1;
end,
end
if ins==0
v=i;
if k>label(v)
k=label(v); v1=v;
end,
end,
end
s(length(s)+1)=v1;
u=v1;
end
min=label(terminal); path(1)=terminal;
i=1;
while path(i)~=start
path(i+1)=f(path(i));
i=i+1 ;
end
path(i)=start;
L=length(path);
path=path(L:-1:1);
测试数据：
[0 2 8 1 Inf Inf Inf Inf Inf Inf Inf;
2 0 6 Inf 1 Inf Inf Inf Inf Inf Inf;
8 6 0 7 5 1 2 Inf Inf Inf Inf;
1 Inf 7 0 Inf Inf 9 Inf Inf Inf Inf;
Inf 1 5 Inf 0 3 Inf 2 9 Inf Inf;
Inf Inf 1 Inf 3 0 4 Inf 6 Inf Inf;
Inf Inf 2 9 Inf 4 0 Inf 3 1 Inf;
Inf Inf Inf Inf 2 Inf Inf 0 7 Inf 9;
Inf Inf Inf Inf 9 6 3 7 0 1 2;
Inf Inf Inf Inf Inf Inf 1 Inf 1 0 4;
Inf Inf Inf Inf Inf Inf Inf 9 2 4 0;]
[0 8 Inf Inf Inf Inf 7 8 Inf Inf Inf;
Inf 0 3 Inf Inf Inf Inf Inf Inf Inf Inf;
Inf Inf 0 5 6 Inf 5 Inf Inf Inf Inf;
Inf Inf Inf 0 1 Inf Inf Inf Inf Inf 12;
Inf Inf 6 Inf 0 2 Inf Inf Inf Inf 10;
Inf Inf Inf Inf 2 0 9 Inf 3 Inf Inf;
Inf Inf Inf Inf Inf 9 0 Inf Inf Inf Inf;
8 Inf Inf Inf Inf Inf Inf 0 9 Inf Inf;
Inf Inf Inf Inf 7 Inf Inf 9 0 2 Inf;
Inf Inf Inf Inf Inf Inf Inf Inf 2 0 2;
Inf Inf Inf Inf 10 Inf Inf Inf Inf Inf 0;];
5.图论_Floryd算法
代码实现：
tulun2.m
a= [ 0,50,inf,40,25,10;
50,0,15,20,inf,25;
inf,15,0,10,20,inf;
40,20,10,0,10,25;
25,inf,20,10,0,55;
10,25,inf,25,55,0];
[D, path]=floyd(a)
floyd.m
function [D,path,min1,path1]=floyd(a,start,terminal)
D=a;n=size(D,1);path=zeros(n,n);
for i=1:n
for j=1:n
if D(i,j)~=inf
path(i,j)=j;
end,
end,
end
for k=1:n
for i=1:n
for j=1:n
if D(i,k)+D(k,j)<D(i,j)
D(i,j)=D(i,k)+D(k,j);
path(i,j)=path(i,k);
end,
end,
end,
end
if nargin==3
min1=D(start,terminal);
m(1)=start;
i=1;
path1=[ ];
while path(m(i),terminal)~=terminal
k=i+1;
m(k)=path(m(i),terminal);
i=i+1;
end
m(i+1)=terminal;
path1=m;
end
测试数据：
[0 8 Inf Inf Inf Inf 7 8 Inf Inf Inf;
Inf 0 3 Inf Inf Inf Inf Inf Inf Inf Inf;
Inf Inf 0 5 6 Inf 5 Inf Inf Inf Inf;
Inf Inf Inf 0 1 Inf Inf Inf Inf Inf 12;
Inf Inf 6 Inf 0 2 Inf Inf Inf Inf 10;
Inf Inf Inf Inf 2 0 9 Inf 3 Inf Inf;
Inf Inf Inf Inf Inf 9 0 Inf Inf Inf Inf;
8 Inf Inf Inf Inf Inf Inf 0 9 Inf Inf;
Inf Inf Inf Inf 7 Inf Inf 9 0 2 Inf;
Inf Inf Inf Inf Inf Inf Inf Inf 2 0 2;
Inf Inf Inf Inf 10 Inf Inf Inf Inf Inf 0;];
6.粒子群_退火：
代码实现：
swap.m
function [ newpath , position ] = swap( oldpath , number )
% 对 oldpath 进 行 互 换 操 作
% number 为 产 生 的 新 路 径 的 个 数
% position 为 对 应 newpath 互 换 的 位 置
m = length( oldpath ) ; % 城 市 的 个 数
newpath = zeros( number , m ) ;
position = sort( randi( m , number , 2 ) , 2 ); % 随 机 产 生 交 换 的 位 置
for i = 1 : number
newpath( i , : ) = oldpath ;
% 交 换 路 径 中 选 中 的 城 市
newpath( i , position( i , 1 ) ) = oldpath( position( i , 2 ) ) ;
newpath( i , position( i , 2 ) ) = oldpath( position( i , 1 ) ) ;
end
pathfare.m
function [ objval ] = pathfare( fare , path )
% 计 算 路 径 path 的 代 价 objval
% path 为 1 到 n 的 排 列 ，代 表 城 市 的 访 问 顺 序 ；
% fare 为 代 价 矩 阵 ， 且 为 方 阵 。
[ m , n ] = size( path ) ;
objval = zeros( 1 , m ) ;
for i = 1 : m
for j = 2 : n
objval( i ) = objval( i ) + fare( path( i , j - 1 ) , path( i , j ) ) ;
end
objval( i ) = objval( i ) + fare( path( i , n ) , path( i , 1 ) ) ;
end
distance.m
function [ fare ] = distance( coord )
% 根 据 各 城 市 的 距 离 坐 标 求 相 互 之 间 的 距 离
% fare 为 各 城 市 的 距 离 ， coord 为 各 城 市 的 坐 标
[ v , m ] = size( coord ) ; % m 为 城 市 的 个 数
fare = zeros( m ) ;
for i = 1 : m % 外 层 为 行
for j = i : m % 内 层 为 列
fare( i , j ) = ( sum( ( coord( : , i ) - coord( : , j ) ) .^ 2 ) ) ^ 0.5 ;
fare( j , i ) = fare( i , j ) ; % 距 离 矩 阵 对 称
end
end
myplot.m
function [ ] = myplot( path , coord , pathfar )
% 做 出 路 径 的 图 形
% path 为 要 做 图 的 路 径 ，coord 为 各 个 城 市 的 坐 标
% pathfar 为 路 径 path 对 应 的 费 用
len = length( path ) ;
clf ;
hold on ;
title( [ '近似最短路径如下，路程为' , num2str( pathfar ) ] ) ;
plot( coord( 1 , : ) , coord( 2 , : ) , 'ok');
pause( 0.4 ) ;
for ii = 2 : len
plot( coord( 1 , path( [ ii - 1 , ii ] ) ) , coord( 2 , path( [ ii - 1 , ii ] ) ) , '-b');
x = sum( coord( 1 , path( [ ii - 1 , ii ] ) ) ) / 2 ;
y = sum( coord( 2 , path( [ ii - 1 , ii ] ) ) ) / 2 ;
text( x , y , [ '(' , num2str( ii - 1 ) , ')' ] ) ;
pause( 0.4 ) ;
end
plot( coord( 1 , path( [ 1 , len ] ) ) , coord( 2 , path( [ 1 , len ] ) ) , '-b' ) ;
x = sum( coord( 1 , path( [ 1 , len ] ) ) ) / 2 ;
y = sum( coord( 2 , path( [ 1 , len ] ) ) ) / 2 ;
text( x , y , [ '(' , num2str( len ) , ')' ] ) ;
pause( 0.4 ) ;
hold off ;
clear;
% 程 序 参 数 设 定
Coord = ... % 城 市 的 坐 标 Coordinates
[ 0.6683 0.6195 0.4 0.2439 0.1707 0.2293 0.5171 0.8732 0.6878 0.8488 ; ...
0.2536 0.2634 0.4439 0.1463 0.2293 0.761 0.9414 0.6536 0.5219 0.3609 ] ;
t0 = 1 ; % 初 温 t0
iLk = 20 ; % 内 循 环 最 大 迭 代 次 数 iLk
oLk = 50 ; % 外 循 环 最 大 迭 代 次 数 oLk
lam = 0.95 ; % λ lambda
istd = 0.001 ; % 若 内 循 环 函 数 值 方 差 小 于 istd 则 停 止
ostd = 0.001 ; % 若 外 循 环 函 数 值 方 差 小 于 ostd 则 停 止
ilen = 5 ; % 内 循 环 保 存 的 目 标 函 数 值 个 数
olen = 5 ; % 外 循 环 保 存 的 目 标 函 数 值 个 数
% 程 序 主 体
m = length( Coord ) ; % 城 市 的 个 数 m
fare = distance( Coord ) ; % 路 径 费 用 fare
path = 1 : m ; % 初 始 路 径 path
pathfar = pathfare( fare , path ) ; % 路 径 费 用 path fare
ores = zeros( 1 , olen ) ; % 外 循 环 保 存 的 目 标 函 数 值
e0 = pathfar ; % 能 量 初 值 e0
t = t0 ; % 温 度 t
for out = 1 : oLk % 外 循 环 模 拟 退 火 过 程
ires = zeros( 1 , ilen ) ; % 内 循 环 保 存 的 目 标 函 数 值
for in = 1 : iLk % 内 循 环 模 拟 热 平 衡 过 程
[ newpath , v ] = swap( path , 1 ) ; % 产 生 新 状 态
e1 = pathfare( fare , newpath ) ; % 新 状 态 能 量
% Metropolis 抽 样 稳 定 准 则
r = min( 1 , exp( - ( e1 - e0 ) / t ) ) ;
if rand < r
path = newpath ; % 更 新 最 佳 状 态
e0 = e1 ;
end
ires = [ ires( 2 : end ) e0 ] ; % 保 存 新 状 态 能 量
% 内 循 环 终 止 准 则 ：连 续 ilen 个 状 态 能 量 波 动 小 于 istd
if std( ires , 1 ) < istd
break ;
end
end
ores = [ ores( 2 : end ) e0 ] ; % 保 存 新 状 态 能 量
% 外 循 环 终 止 准 则 ：连 续 olen 个 状 态 能 量 波 动 小 于 ostd
if std( ores , 1 ) < ostd
break ;
end
t = lam * t ;
end
pathfar = e0 ;
% 输 入 结 果
fprintf( '近似最优路径为：\n ' )
%disp( char( [ path , path(1) ] + 64 ) ) ;
disp(path)
fprintf( '近似最优路径路程\tpathfare=' ) ;
disp( pathfar ) ;
myplot( path , Coord , pathfar ) ;测试数据：
Coord = [ 66.83 61.95 40 24.39 17.07 22.93 51.71 87.32 68.78 84.88 50 40 25 ;
25.36 26.34 44.39 14.63 22.93 76.1 94.14 65.36 52.19 36.09 30 20 26] ;7.粒子群_种群竞争
代码实现：
fun.m:
function dx=fun(t,x,r1,r2,n1,n2,s1,s2)
r1=1;
r2=1;
n1=100;
n2=100;
s1=0.5;
s2=2;
dx=[r1*x(1)*(1-x(1)/n1-s1*x(2)/n2);r2*x(2)*(1-s2*x(1)/n1-x(2)/n2)];
p3.m:
h=0.1;%所取时间点间隔
ts=[0:h:30];%时间区间
x0=[10,10];%初始条件
opt=odeset('reltol',1e-6,'abstol',1e-9);%相对误差1e-6，绝对误差1e-9
[t,x]=ode45(@fun,ts,x0,opt);%使用5级4阶龙格—库塔公式计算
plot(t,x(:,1),'r',t,x(:,2),'b','LineWidth',2),grid;
pause;
plot(x(:,1),x(:,2),'LineWidth',2),grid %作相轨线
8.排队论：
代码1：
clear
clc
%*****************************************
%初始化顾客源
%*****************************************
%总仿真时间
Total_time = 10;
%队列最大长度
N = 10000000000;
%到达率与服务率
lambda = 10;
mu = 6;
%平均到达时间与平均服务时间
arr_mean = 1/lambda;
ser_mean = 1/mu;
arr_num = round(Total_time*lambda*2);
events = [];
%按负指数分布产生各顾客达到时间间隔
events(1,:) = exprnd(arr_mean,1,arr_num);
%各顾客的到达时刻等于时间间隔的累积和
events(1,:) = cumsum(events(1,:));
%按负指数分布产生各顾客服务时间
events(2,:) = exprnd(ser_mean,1,arr_num);
%计算仿真顾客个数，即到达时刻在仿真时间内的顾客数
len_sim = sum(events(1,:)<= Total_time);
%*****************************************
%计算第 1个顾客的信息
%*****************************************
%第 1个顾客进入系统后直接接受服务，无需等待
events(3,1) = 0;
%其离开时刻等于其到达时刻与服务时间之和
events(4,1) = events(1,1)+events(2,1);
%其肯定被系统接纳，此时系统内共有
%1个顾客，故标志位置1
events(5,1) = 1;
%其进入系统后，系统内已有成员序号为 1
member = [1];
for i = 2:arr_num
%如果第 i个顾客的到达时间超过了仿真时间，则跳出循环
if events(1,i)>Total_time
break;
else
number = sum(events(4,member) > events(1,i));
%如果系统已满，则系统拒绝第 i个顾客，其标志位置 0
if number >= N+1
events(5,i) = 0;
%如果系统为空，则第 i个顾客直接接受服务
else
if number == 0
%其等待时间为 0
2009.1516
%PROGRAMLANGUAGEPROGRAMLANGUAGE
events(3,i) = 0;
%其离开时刻等于到达时刻与服务时间之和
events(4,i) = events(1,i)+events(2,i);
%其标志位置 1
events(5,i) = 1;
member = [member,i];
%如果系统有顾客正在接受服务，且系统等待队列未满，则 第 i个顾客进入系统
else len_mem = length(member);
%其等待时间等于队列中前一个顾客的离开时刻减去其到 达时刻
events(3,i)=events(4,member(len_mem))-events(1,i);
%其离开时刻等于队列中前一个顾客的离开时刻加上其服
%务时间
events(4,i)=events(4,member(len_mem))+events(2,i);
%标识位表示其进入系统后，系统内共有的顾客数
events(5,i) = number+1;
member = [member,i];
end
end
end
end
%仿真结束时，进入系统的总顾客数
len_mem = length(member);
%*****************************************
%输出结果
%*****************************************
%绘制在仿真时间内，进入系统的所有顾客的到达时刻和离
%开时刻曲线图（stairs：绘制二维阶梯图）
stairs([0 events(1,member)],0:len_mem);
hold on;
stairs([0 events(4,member)],0:len_mem,'.-r');
legend('到达时间 ','离开时间 ');
hold off;
grid on;
%绘制在仿真时间内，进入系统的所有顾客的停留时间和等
%待时间曲线图（plot：绘制二维线性图）
figure;
plot(1:len_mem,events(3,member),'r-*',1: len_mem,events(2,member)+events(3,member),'k-');
legend('等待时间 ','停留时间 ');
grid on;代码2：
s=2;
mu=4;
lambda=3;
ro=lambda/mu;
ros=ro/s;
sum1=0;
for i=0:(s-1)
sum1=sum1+ro.^i/factorial(i);
end
sum2=ro.^s/factorial(s)/(1-ros);
p0=1/(sum1+sum2);
p=ro.^s.*p0/factorial(s)/(1-ros);
Lq=p.*ros/(1-ros);
L=Lq+ro;
W=L/lambda;
Wq=Lq/lambda;
fprintf('排队等待的平均人数为%5.2f人\n',Lq)
fprintf('系统内的平均人数为%5.2f人\n',L)
fprintf('平均逗留时间为%5.2f分钟\n',W*60)
fprintf('平均等待时间为%5.2f分种\n',Wq*60)
扫一扫
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
