matlab:数学模型总结（一）

xiaoxiao2023-10-26 152

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.多属性决策：

代码实现：

测试数据：

[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)

最新回复(0)