数学建模常用算法Matlab代码_数学建模常用的matlab代码

目录

1.通用Topsis算法

2.通用Floyd算法

3.通用Dijkstra代码

4.层次分析法代码

5.插值与拟合代码

6.灰色关联分析代码

7.主成分分析代码

8.蒙特卡洛算法代码

9.规划模型代码

10.灰色预测算法代码

11.智能算法之模拟退火算法

---

1.通用Topsis算法

%西南石油计科罗文睿
[n,m]=size(X);
disp(['共有' num2str(n) '个评价对象,' num2str(m) '个评价指标'])
Judge=input(['这' num2str(m) '个评价指标需要经过正向化处理，需要请输入1，不需要请输入0：'])
if Judge ==1
    Position=input('请输入需要正向化处理的指标所在的列，如2，3，6列，[2,3,6]：');
    disp('请输入需要处理的这些列的指标类型：1.极小型 2.中间型 3.区间型')
    Type = input('2列是极小型，3列是中间型，6列是区间型,[2,3,6]：');
 
    for i=1:size(Position,2)
        X(:,Position(i))=Positivization(X(:,Position(i)),Type(i),Position(i));
    end
    disp('正向化后的矩阵X= ')
    disp(X)
end
 
Z=X./repmat(sum(X.*X).^0.5,n,1);
disp('标准化矩阵Z= ')
disp(Z)
 
D_P=sum((Z-repmat(max(Z),n,1)).^2,2).^0.5;
D_N=sum((Z-repmat(min(Z),n,1)).^2,2).^0.5;
S=D_N./(D_P+D_N);
disp('最后的得分为：')
stand_S=S/sum(S)
[sorted_S,index]=sort(stand_S,'descend')
 
%以下为自定义函数，需单独存放在一个函数里
function [posit_x]=Positivization(x,type,i)
    if type==1
        disp(['第' num2str(i) '列是极小型，正在正向化'])
        posit_x=Min2Max(x);
        disp(['第' num2str(i) '列极小型正向化处理完成'])
        disp('-------------------分界线----------------------')
    elseif type==2
        disp(['第' num2str(i) '列是中间型'])
        best=input('请输入最佳的那一个值：')
        posit_x=Mid2Max(x,best);
        disp(['第' num2str(i) '列中间型正向化处理完成'])
        disp('------------------分界线-----------------------')
    elseif type==3
        disp(['第' num2str(i) '列是区间型'])
        a=input('请输入区间的下界：');
        b=input('请输入区间的上界：');
        posit_x=lnter2Max(x,a,b);
        disp(['第' num2str(i) '列区间型正向化处理完成'])
        disp('------------------分界线-----------------------')
    else
        disp('没有这种类型的指标，请检查Type向量中是否有除了1，2，3之外的其他值')
    end
end
 
%自定义函数 极小型转极大型
function [posit_x]=Min2Max(x)
    posit_x=max(x)-x;
end
 
%自定义函数 中间型转极大型
function [posit_x]=Mid2Max(x,best)
    M=max(abs(x-best));
    posit_x=1-abs(x-best)/M;
end
 
%自定义函数 区间型转极大型
function [posit_x]=lnter2Max(x,a,b)
    r_x=size(x,1);
    M=max([a-min(x),max(x)-b]);
    posit_x=zeros(r_x,1);
    for i=1:r_x
        if x(i)<a
            posit_x(i)=1-(a-x(i))/M;
        elseif x(i)>b
            posit_x(i)=1-(x(i)-b)/M;
        else
            posit_x(i)=1;
        end
    end
end

2.通用Floyd算法

%西南石油计科罗文睿
function [dist,mypath]=myfloyd(a,sb,db);
%a-邻接矩阵 元素a(i,j)-顶点i到j之间的直达距离，可以是有向的
%sb-起点的标号 db-终点的标号
%输出：dist-最短路的距离 mypath-最短路的路径
n=size(a,1);
path=zeros(n);
for k=1:n
    for i=1:n
        for j=1:n
            if a(i,j)>a(i,k)+a(k,j)
                a(i,j)=a(i,k)+a(k,j);
                path(i,j)=k;
            end
        end
    end
end
dist=a(sb,db);
parent=path(sb,:);
parent(parent==0)=sb;
mypath=db;
t=db;
while t~=sb
    p=parent(t);
    mypath=[p,mypath];
    t=p;
end

3.通用Dijkstra代码

%西南石油计科罗文睿
function [mydistance,mypath]=mydijkstra(a,sb,db);
%a是邻接矩阵 a(i,j)-i到j之间的距离，可以是有向的
%sb-起点的标号 db-终点的标号
%mydistance-最短路的距离 mypath-最短路的路径
n=size(a,1);
visited(1:n)=0;
distance(1:n)=inf;
distance(sb)=0;
visited(sb)=1;
u=sb;
parent(1:n)=0;
for i=1:n-1
    id=find(visited==0);
    for v=id
        if a(u,v)+distance(u)<distance(v)
            distance(v)=distance(u)+a(u,v);
            parent(v)=u;
        end
    end
    temp=distance;
    temp(visited==1)=inf;
    [t,u]=min(temp);
    visited(u)=1;
end
mypath=[];
if parent(db)~=0
    t=db;
    mypath=[db];
    while t~=sb
        P=parent(t);
        mypath=[P mypath];
        t=P;
    end
end
mydistance=distance(db);

4.层次分析法代码

%西南石油计科罗文睿
% 使用方法
%（1）构造判断矩阵A
%（2）将下文代码复制粘贴到Matlab中即可
% 例如：A=[1 3 5;0.33 1 3;0.2 0.33,1]
 
disp('请输入准则层判断矩阵A(n阶)');
A=input('A=');
[n,n]=size(A);
[V,D]=eig(A);%求得特征向量和特征值
            %求出最大特征值和它所对应的特征向量
tempNum=D(1,1);
pos=1;
for h=1:n
    if D(h,h)>tempNum
        tempNum=D(h,h);
        pos=h;
    end
end    
w=abs(V(:,pos));
w=w/sum(w);
t=D(pos,pos);
disp('准则层特征向量w=');disp(w);disp('准则层最大特征根t=');disp(t);
         %以下是一致性检验
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 1.60 1.61 1.615 1.62 1.63];
CR=CI/RI(n);
if CR<0.10
    disp('此矩阵的一致性可以接受!');
    disp('CI=');disp(CI);
    disp('CR=');disp(CR);
else disp('此矩阵的一致性验证失败，请重新进行评分!');
end
 
 
disp('请输入方案层各因素对准则层各因素权重的成对比较阵');
for i=1:n
    disp('请输入第');disp(i);disp('个准则层因素的判断矩阵B');disp(i);

5.插值与拟合代码

% 一维插值步骤
%（1）输入已知数据，x,y
%（2）输入待插自变量的值x1
x=1:12;
y=[5 8 9 15 25 29 31 30 22 25 27 24];
x1=1:0.1:12;
t=interp1(x,y,x1,'spline');%
 
plot(x1,t,'r:')   %作图
 
xlabel('x'),ylabel('y')
 
 
% 二维插值步骤
%（1）先输入二维数据的x,y坐标值
%（2）输入Z数据
%（3）输入待插点的x,y坐标
%（4）应用函数插值即可
x=1:5;
 
y=1:3;
 
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];
 
mesh(x,y,temps);
xi=1:0.2:5;
 
yi=1:0.2:3;
 
zi=interp2(x,y,temps,xi',yi,'cubic');
 
mesh(xi,yi,zi);
 
% 多项式拟合步骤
% (1)输入待拟合数据x,y
% (2)输入函数公式进行拟合
x=0:0.1:1;
 
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66   9.56 9.48 9.30 11.2];
 
A=polyfit(x,y,2)
 
% 指定函数拟合步骤
% (1)输入待拟合数据x,y
%（2）指定函数关系式
syms t;
x=[0;0.4;1.2;2;2.8;3.6;4.4;5.2;6;7.2;8;9.2;10.4;11.6;12.4;13.6;14.4;15];
y=[1;0.85;0.29;-0.27;-0.53;-0.4;-0.12;0.17;0.28;0.15;-0.03;-0.15;-0.071;0.059;0.09;0.032;-0.015;-0.02];%指定函数形式为f(t)=acos(kt)e^(wt),进行拟合
f=fittype('a*cos(k*t)*exp(w*t)','independent','t','coefficients',{'a','k','w'});
cfun=fit(x,y,f)     %显示拟合函数
xi=0:.1:20;
yi=cfun(xi);
plot(x,y,'r*',xi,yi,'b-');
 

6.灰色关联分析代码

%灰色关联分析步骤
%【1】确定比较对象（评价对象）(就是数据，并且需要进行规范化处理，就是标准化处理，见下面例题的表格数据)和参考数列（评价标准，一般该列数列都是1，就是最优的的情况） 
%【2】确定各个指标权重，可用层次分析确定 
%【3】计算灰色关联系数 
%【4】计算灰色加权关联度 
%【5】评价分析
 
x1=[1.14 1.49 1.69 2.12 2.43 4.32 5.92 6.07 7.85;3.30 3.47 3.61 3.80 4.00 4.19 4.42 4.61 4.80;6.00 6.00 6.00 7.50 7.50 7.50 9.00 9.00 9.00;1.20 1.20 1.80 1.80 1.80 2.40 2.70 3.60 4.00;4.87 5.89 6.76 7.97 8.84 10.05 11.31 12.25 11.64]%原始数据5行9列
x=x1;
for i=1:5
    for j=1:9
  x(i,j)=x(i,j)/x1(1,j)
 end
end
x1=x
for i=1:5
    for j=1:9
  x(i,j)=abs(x(i,j)-x1(i,1))
 end
end
max=x(1,1)
min=x(1,1)
for i=1:5
    for j=1:9
 if x(i,j)>=max
     max=x(i,j)
 end
  end
end
for i=1:5
    for j=1:9
 if x(i,j)<=min
     min=x(i,j)
  end
  end
end
k=0.5 %分辨系数取值
l=(min+k*max)./(x+k*max)%求关联系数矩阵
guanliandu=sum(l')/n
[rs,rind]=sort(guanliandu,'descend') %对关联度进行排序

7.主成分分析代码

% PCA步骤：
 
%（1）对原始数据进行标准化处理
 
%（2）计算样本相关系数矩阵
 
%（3）计算相关系数矩阵R的特征值和相应的特征向量
 
%（4）选择重要的主成分，写出主成分表达式
 
%下例中企业综合实力排序问题，其中各列分别为：企业序号；净利润率；固定资产利润率；总产值利润率；销售收入利润率；产品成本利润率；物耗利润率；人均利润；流动资金
 
x =
 
    1.0000   40.4000   24.7000    7.2000    6.1000    8.3000    8.7000    2.4420   20.0000
    2.0000   25.0000   12.7000   11.2000   11.0000   12.9000   20.2000    3.5420    9.1000
    3.0000   13.2000    3.3000    3.9000    4.3000    4.4000    5.5000    0.5780    3.6000
    4.0000   22.3000    6.7000    5.6000    3.7000    6.0000    7.4000    0.1760    7.3000
    5.0000   34.3000   11.8000    7.1000    7.1000    8.0000    8.9000    1.7260   27.5000
    6.0000   35.6000   12.5000   16.4000   16.7000   22.8000   29.3000    3.0170   26.6000
    7.0000   22.0000    7.8000    9.9000   10.2000   12.6000   17.6000    0.8470   10.6000
    8.0000   48.4000   13.4000   10.9000    9.9000   10.9000   13.9000    1.7720   17.8000
    9.0000   40.6000   19.1000   19.8000   19.0000   29.7000   39.6000    2.4490   35.8000
   10.0000   24.8000    8.0000    9.8000    8.9000   11.9000   16.2000    0.7890   13.7000
   11.0000   12.5000    9.7000    4.2000    4.2000    4.6000    6.5000    0.8740    3.9000
   12.0000    1.8000    0.6000    0.7000    0.7000    0.8000    1.1000    0.0560    1.0000
   13.0000   32.3000   13.9000    9.4000    8.3000    9.8000   13.3000    2.1260   17.1000
   14.0000   38.5000    9.1000   11.3000    9.5000   12.2000   16.4000    1.3270   11.6000
   15.0000   26.2000   10.1000    5.6000   15.6000    7.7000   30.1000    0.1260   25.9000
 
A=x;
a=size(A,1);%获得矩阵A的行大小
b=size(A,2);%获得矩阵A的列大小
for i=1:b
    SA(:,i)=(A(:,i)-mean(A(:,i)))/std(A(:,i));%std函数是用来求向量的标准差
end
% %计算相关系数矩阵的特征值和特征向量
CM=corrcoef(SA);%计算相关系数矩阵
[V,D]=eig(CM);%计算特征值和特征向量
for j=1:b
    DS(j,1)=D(b+1-j,b+1-j);%对特征值按降序排列
end
for i=1:b
    DS(i,2)=DS(i,1)/sum(DS(:,1));%贡献率
    DS(i,3)=sum(DS(1:i,1))/sum(DS(:,1));%累计贡献率
end
% % 选择主成分及对应的特征向量
T=0.9;%主成分信息保留率
for k=1:b
    if DS(k,3)>=T
        Com_num=k;
        break;
    end
end
%提取主成分对应的特征向量
for j=1:Com_num
    PV(:,j)=V(:,b+1-j);
end
% % 计算各评价对象的主成分得分
new_score=SA*PV;
for i=1:a
    total_score(i,1)=sum(new_score(i,:));
    total_score(i,2)=i;
end
result_report=[new_score,total_score];%将各主成分得分与总分放在同一个矩阵中
result_report=sortrows(result_report,-4);%按总分降序排序
% % 输出模型及结果报告
disp('特征值及其贡献率，累加贡献率：')
DS
disp('信息保留率T对应的主成分数与特征向量：')
Com_num
PV
disp('主成分得分及排序（按第四列的总分进行排序，前三列为个主成分得分，第五列为企业编号）')
result_report

8.蒙特卡洛算法代码

%蒙特卡洛法是经过大量事件的统计结果来实现一些确定性问题的计算。
%使用蒙特卡洛法必须使用计算机生成相关分布的随机数。
%例如：y = x^2 ，y = 12 - x与X轴在第一象限与X轴围成一个曲边三角形。设计一个随机试验，求该图形的近似值。
x=0:0.25:12
y1=x.^2;
y2=12-x;
plot(x,y1,x,y2)
xlabel('x');ylabel('y');
legend('y1=x^2','y2=12-x');
title('罗文睿绘制');
axis([0 15 0 15]);
text(3,9,'交点');
grid on
 
x=unifrnd(0,12,[1,10000000]);
y=unifrnd(0,9,[1,10000000]);
frequency=sum(y<x.^2 & x<=3)+sum(y<12-x & x>=3)
area=12*9*frequency/10^7

9.规划模型代码

%线性规划步骤
%（1）定义决策变量
%（1）构造目标函数
%（2）寻找限制条件
%（4）按照步骤带入Matlab函数
 
c=-[5,3]';  
   
A=[2,1;1,2];  
   
b=[40,50]'; 
   
L=[0, 0];
   
[x,fmin]=linprog(c,A,b,[],[],L);
   
Pmax=-fmin
   x1=x(1),  x2=x(2)
 
 
 
% 非线性规划代码
H=[1 -1; -1 2]; 
      c=[-2 ;-6];A=[1 1; -1 2];b=[2;2];
      Aeq=[];beq=[]; VLB=[0;0];VUB=[];
      [x,z]=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)
 
% 多目标规划代码
 
f=[3;-2]; 
A=[2,3;2,1]; 
b=[18;10]; 
lb=[0;0];
[x,fval]=linprog(f,A,b,[],[],lb)
 
f=[-4;-3]; 
A=[2,3;2,1]; 
b=[18;10]; 
lb=[0;0];
[x,fval]=linprog(f,A,b,[],[],lb)

10.灰色预测算法代码

% 灰色预测步骤
%（1）输入前期的小样本数据
%（2）输入预测个数
%（3）运行
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;
t_test=input('输入需要预测的个数');
i=1:t_test+n;
yys(i+1)=(y(1)-t).*exp(-a.*i)+t;
yys(1)=y(1);
for j=n+t_test:-1:2
    ys(j)=yys(j)-yys(j-1);
end
x=1:n;
xs=2:n+t_test;
yn=ys(2:n+t_test);
plot(x,y,'^r',xs,yn,'*-b');
det=0;
for i=2:n
    det=det+abs(yn(i)-y(i));
end
det=det/(n-1);
disp(['百分绝对误差为：',num2str(det),'%']);
disp(['预测值为：',num2str(ys(n+1:n+t_test))]);

11.智能算法之模拟退火算法

 
%生成初始解，求目标函数f(x)=x1^2+x2^2+8在x1^2-x2>0;-x1-x2^2+2=0约束下的最小值问题  
sol_new2=1;%（1）解空间（初始解）  
sol_new1=2-sol_new2^2;  
sol_current1 = sol_new1;   
sol_best1 = sol_new1;  
sol_current2 = sol_new2;   
sol_best2 = sol_new2;  
E_current = inf;  
E_best = inf;  
  
rand('state',sum(clock)); %初始化随机数发生器  
t=90; %初始温度  
tf=89.9; %结束温度  
a = 0.99; %温度下降比例  
  
while t>=tf%（7）结束条件  
    for r=1:1000 %退火次数  
          
        %产生随机扰动（3）新解的产生  
        sol_new2=sol_new2+rand*0.2;  
        sol_new1=2-sol_new2^2;  
          
        %检查是否满足约束  
        if sol_new1^2-sol_new2>=0 && -sol_new1-sol_new2^2+2==0 && sol_new1>=0 &&sol_new2>=0  
        else  
            sol_new2=rand*2;  
            sol_new1=2-sol_new2^2;  
            continue;  
        end  
          
        %退火过程  
        E_new=sol_new1^2+sol_new2^2+8;%（2）目标函数  
        if E_new<E_current%（5）接受准则  
                E_current=E_new;  
                sol_current1=sol_new1;  
                sol_current2=sol_new2;  
                if E_new<E_best  
                    %把冷却过程中最好的解保存下来  
                    E_best=E_new;  
                    sol_best1=sol_new1;  
                    sol_best2=sol_new2;  
                end  
        else  
                if rand<exp(-(E_new-E_current)/t)%（4）代价函数差  
                    E_current=E_new;  
                    sol_current1=sol_new1;  
                    sol_current2=sol_new2;  
                else  
                    sol_new1=sol_current1;  
                    sol_new2=sol_current2;  
                end  
        end  
        plot(r,E_best,'*')  
        hold on  
    end  
    t=t*a;%（6）降温  
end  
  
disp('最优解为：')  
disp(sol_best1)  
disp(sol_best2)  
disp('目标表达式的最小值等于：')  
disp(E_best)  

未完待续