目录
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)
未完待续