%原始数据标准化处理
clc;
clear all;
load x1.txt;
load x2.txt;
load v1.txt;
[n,p]=size(x1);
[n1,p1]=size(v1);
z=zeros(n,p); %标准化矩阵z
x_ave=zeros(1,p); %原始数据矩阵x的各个指标的平均值x_zve
s=zeros(1,p); %原始数据矩阵x的各个指标的标准差s
v=zeros(p,p); %相关性矩阵r
z_ave=zeros(1,p); %标准化矩阵z的各个指标的平均值z_ave
%求原始数据各个指标的平均值x_ave
for j=1:p
for i=1:n
x_ave(1,j)=x_ave(1,j)+x1(i,j);
end
x_ave(1,j)=x_ave(1,j)/n;
end
%计算原始数据矩阵的各个指标的标准差s
for j=1:p
for i=1:n
s(1,j)=s(1,j)+(x1(i,j)-x_ave(1,j))^2;
end
s(1,j)=sqrt(s(1,j)/n);
end
%求标准化处理后的数据矩阵z
for i=1:n
for j=1:p
z(i,j)=(x1(i,j)-x_ave(1,j))/s(1,j);
end
end
%z=abs(z);
%求标准化后的数据矩阵中各个指标的平均值z_ave
for j=1:p
for i=1:n
z_ave(1,j)=z_ave(1,j)+z(i,j);
end
z_ave(1,j)=z_ave(1,j)/n;
end
%计算相关性矩阵v
a=0;
b=0;
c=0;
for k=1:p
for j=1:p
for i=1:n
a=a+(z(i,k)-z_ave(1,k))*(z(i,j)-z_ave(1,j));
b=b+(z(i,k)-z_ave(1,k))^2;
c=c+(z(i,j)-z_ave(1,j))^2;
end
b=sqrt(b);
c=sqrt(c);
v(k,j)=a/(b*c);
a=0;
b=0;
c=0;
end
end
%求相关性矩阵的特征值矩阵和特征向量矩阵
[V,D]=eig(v); % V是特征向量矩阵,D是特征值矩阵
TD=diag(D); % TD是特征值组成的列向量
%确定总的方差贡献率
rp=0;
for i=1:p
rp=rp+TD(i,1);
end
%将特征值从大到小排序
%TDF=zeros(p,1);
%VF=zeros(p,p);
%for i=1:p
% TDF(p+1-i,1)=TD(i,1);
%end
%for i=1:p
% for j=1:p
% VF(j,p+1-i)=V(j,i);
% end
%end
[TDF,I]=sort(TD,'descend'); %将特征值TD由大到小排序
VF=V(:,I); %特征向量与特征值相对应
%确定各个主成分贡献率(即各个主成分的权重)
w=zeros(p,1);
for i=1:p
w(i,1)=TDF(i,1)/rp;
end
%确定m个主元
m=0;
for k=1:p
rt=0;
for t=1:k
rt=rt+w(t,1);
end
if rt<=0.90
m=m+1;
end
end
%确定m个主元分别对应的p个特征向量矩阵T
%T=zeros(p,m);
%TF=zeros(p,m);
%for i=1:p
% for j=1:m
% TF(i,j)=VF(i,j);
% end
%end
%for i=1:m
% for j=1:p
% T(j,m+1-i)=TF(j,i);
% end
%end
%求主元综合样本矩阵
%F=zeros(n,m);
%F=z*T;
%求主元综合样本矩阵F
for i=1:m
X1(:,i)=z*VF(:,i);
end
F=0;
for i=1:m
F=F+X1(:,i)*VF(:,i)';
end
%X1=z*VF;
%建立多项指标综合评价数学模型
f=zeros(n,1);
for i=1:n
for j=1:m
f(i,1)=f(i,1)+w(j,1)*X1(i,j);
end
f(i,1)=f(i,1);
end
f=abs(f);
g=1:25;
%figure(1);
%plot(g,f,'-r.','MarkerSize',20);
%xlabel('测试点');
%ylabel('各个变电站所有测试点的腐蚀分值');
%figure(2);
%plot(g,v1,'-r.','MarkerSize',20);
%xlabel('测试点');
%ylabel('各个变电站所有测试点的腐蚀速率');
%figure(3);
%plot(g,x2,'-r.','MarkerSize',20);
%xlabel('测试点');
%ylabel('各个变电站所有测试点的电导率');
figure(4);
subplot(2,2,1:2);
plot(g(2:5),f(2:5),'-k.','MarkerSize',10); hold on;
plot(g(6:9),f(6:9),'-k.','MarkerSize',10); hold on;
plot(g(10:12),f(10:12),'-k.','MarkerSize',10); hold on;
plot(g(14:15),f(14:15),'-k.','MarkerSize',10); hold on;
plot(g(16:17),f(16:17),'-k.','MarkerSize',10); hold on;
plot(g(19:20),f(19:20),'-k.','MarkerSize',10); hold on;
plot(g(22:23),f(22:23),'-k.','MarkerSize',10); hold on;
plot(g(24:25),f(24:25),'-k.','MarkerSize',10); hold on;
xlabel('a) 测试点');
ylabel('腐蚀分值');
%figure(5);
subplot(2,2,3);
plot(g(2:5),v1(2:5),'-k.','MarkerSize',10); hold on;
plot(g(6:9),v1(6:9),'-k.','MarkerSize',10); hold on;
plot(g(10:12),v1(10:12),'-k.','MarkerSize',10); hold on;
plot(g(14:15),v1(14:15),'-k.','MarkerSize',10); hold on;
plot(g(16:17),v1(16:17),'-k.','MarkerSize',10); hold on;
plot(g(19:20),v1(19:20),'-k.','MarkerSize',10); hold on;
plot(g(22:23),v1(22:23),'-k.','MarkerSize',10); hold on;
plot(g(24:25),v1(24:25),'-k.','MarkerSize',10); hold on;
xlabel('b) 测试点');
ylabel('腐蚀速率');
%figure(6);
subplot(2,2,4);
plot(g(2:5),x2(2:5),'-k.','MarkerSize',10); hold on;
plot(g(6:9),x2(6:9),'-k.','MarkerSize',10); hold on;
plot(g(10:12),x2(10:12),'-k.','MarkerSize',10); hold on;
plot(g(14:15),x2(14:15),'-k.','MarkerSize',10); hold on;
plot(g(16:17),x2(16:17),'-k.','MarkerSize',10); hold on;
plot(g(19:20),x2(19:20),'-k.','MarkerSize',10); hold on;
plot(g(22:23),x2(22:23),'-k.','MarkerSize',10); hold on;
plot(g(24:25),x2(24:25),'-k.','MarkerSize',10); hold on;
xlabel('c) 测试点');
ylabel('电导率');