数学建模——熵权+TOPSIS+肘部法则+系统聚类

一、起因

我本科的数学建模队长找上我，让我帮她写下matlab代码，当然用的模型还是曾经打比赛的模型，所以虽然代码量多，但是写的很快，也是正逢中秋，有点时间。
当然我也没想到，研一刚开学这论文就这么水灵灵的要发出来了，队长的实力还是遥遥领先啊，在这里，祝她们的论文一次见刊！

二、代码展示

差不多有460行代码这样

%% 获取数据
[num,txt,raw] = xlsread('C:\Users\19468\Desktop\zhibiao1.xlsx',1);%% 标准化MIN = min(num); %各指标的最大值
MAX = max(num); %各指标的最小值
NUM_B = [];     %保存标准化后的数据
for i = 1:size(num,1)    % 运行31个省份for j = 1:size(num,2) % 1个省份的23个数据NUM_B(i,j) = (num(i,j)-MIN(j))/(MAX(j) - MIN(j));end
end%% 对个别指标正向化
% 对指标在第5，10，11，13，14、22列进行正向化
Z = [5,10,11,13,14,22]; % 保存正向化列数
NUM_B_Z = NUM_B;     % 存放标准化正向化后的数据
for i = 1:size(num,1)    % 运行31个省份for j=1:size(Z,2)    % 5列数据NUM_B_Z(i,Z(j)) = 1 - NUM_B(i,Z(j));end
end%% 对二级指标进行熵权+TOPSIS，保留各省份各一级指标的数据
TOP_one   = [1,2,3,4,5];       % 第一个一级指标
TOP_two   = [6,7,8];           % 第二个一级指标
TOP_three = [9,10,11,12,13];   % 第三个一级指标
TOP_four  = [14,15];           % 第四个一级指标
TOP_five  = [16,17,18,19];     % 第五个一级指标
TOP_six   = [20,21,22,23];     % 第六个一级指标TOP_one_data   = NUM_B_Z(:,1:5);    % 第一个一级指标的数据
TOP_two_data   = NUM_B_Z(:,6:8);    % 第二个一级指标的数据
TOP_three_data = NUM_B_Z(:,9:13);   % 第三个一级指标的数据
TOP_four_data  = NUM_B_Z(:,14:15);  % 第四个一级指标的数据
TOP_five_data  = NUM_B_Z(:,16:19);  % 第五个一级指标的数据
TOP_six_data   = NUM_B_Z(:,20:23);  % 第六个一级指标的数据% NUM_B_Z_W(:,:,1) = TOP_one_data;
% NUM_B_Z_W(:,:,2) = TOP_two_data;
% NUM_B_Z_W(:,:,3) = TOP_three_data;
% NUM_B_Z_W(:,:,4) = TOP_four_data;
% NUM_B_Z_W(:,:,5) = TOP_five_data;% 熵权法+TOPSIS
% 指标一
Z = TOP_one_data;
[n,m] = size(Z);
D = zeros(1,m);  % 初始化保存信息效用值的行向量
for i = 1:mx = Z(:,i);  % 取出第i列的指标p = x / sum(x);% 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以要判断n = length(p);   % 向量的长度lnp = zeros(n,1);   % 初始化最后的结果for j = 1:n   % 开始循环if p(j) == 0   % 如果第i个元素为0lnp(j) = 0;  % 那么返回的第i个结果也为0elselnp(j) = log(p(j));endend%到这里结束e = -sum(p .* lnp) / log(n); % 计算信息熵D(i) = 1- e; % 计算信息效用值
end
W = D ./ sum(D);  % 将信息效用值归一化，得到权重
WD_min = min(TOP_one_data);
D_max = max(TOP_one_data);zonghedefen_one = 1:31;  %生成记录点
for i = 1:31%利用（C = D-/((D-)+(D+))）topsis模型与熵权进行结合D_jia  = sqrt(W(1)*(TOP_one_data(i,1)-D_max(1))^2 + W(2)*(TOP_one_data(i,2)-D_max(2))^2 + ...W(3)*(TOP_one_data(i,3)-D_max(3))^2 + W(4)*(TOP_one_data(i,4)-D_max(4))^2 + ...W(5)*(TOP_one_data(i,5)-D_max(5))^2);D_jian = sqrt(W(1)*(TOP_one_data(i,1)-D_min(1))^2 + W(2)*(TOP_one_data(i,2)-D_min(2))^2 + ...W(3)*(TOP_one_data(i,3)-D_min(1))^2 + W(4)*(TOP_one_data(i,4)-D_min(2))^2 + ...W(5)*(TOP_one_data(i,5)-D_min(1))^2);zonghedefen_one(i) = D_jian/(D_jia+D_jian);  % 一次就是产生31个数据
endzonghedefen_one = zonghedefen_one';  % 指标二
Z = TOP_two_data;
[n,m] = size(Z);
D = zeros(1,m);  % 初始化保存信息效用值的行向量
for i = 1:mx = Z(:,i);  % 取出第i列的指标p = x / sum(x);% 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以要判断n = length(p);   % 向量的长度lnp = zeros(n,1);   % 初始化最后的结果for j = 1:n   % 开始循环if p(j) == 0   % 如果第i个元素为0lnp(j) = 0;  % 那么返回的第i个结果也为0elselnp(j) = log(p(j));endend%到这里结束e = -sum(p .* lnp) / log(n); % 计算信息熵D(i) = 1- e; % 计算信息效用值
end
W = D ./ sum(D);  % 将信息效用值归一化，得到权重
WD_min = min(TOP_two_data);
D_max = max(TOP_two_data);zonghedefen_two = 1:31;  %生成记录点
for i = 1:31%利用（C = D-/((D-)+(D+))）topsis模型与熵权进行结合D_jia  = sqrt(W(1)*(TOP_two_data(i,1)-D_max(1))^2 + W(2)*(TOP_two_data(i,2)-D_max(2))^2 + ...W(3)*(TOP_two_data(i,3)-D_max(3))^2);D_jian = sqrt(W(1)*(TOP_two_data(i,1)-D_min(1))^2 + W(2)*(TOP_two_data(i,2)-D_min(2))^2 + ...W(3)*(TOP_two_data(i,3)-D_min(1))^2);zonghedefen_two(i) = D_jian/(D_jia+D_jian);  % 一次就是产生31个数据
endzonghedefen_two = zonghedefen_two'; % 指标三
Z = TOP_three_data;
[n,m] = size(Z);
D = zeros(1,m);  % 初始化保存信息效用值的行向量
for i = 1:mx = Z(:,i);  % 取出第i列的指标p = x / sum(x);% 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以要判断n = length(p);   % 向量的长度lnp = zeros(n,1);   % 初始化最后的结果for j = 1:n   % 开始循环if p(j) == 0   % 如果第i个元素为0lnp(j) = 0;  % 那么返回的第i个结果也为0elselnp(j) = log(p(j));endend%到这里结束e = -sum(p .* lnp) / log(n); % 计算信息熵D(i) = 1- e; % 计算信息效用值
end
W = D ./ sum(D);  % 将信息效用值归一化，得到权重
WD_min = min(TOP_three_data);
D_max = max(TOP_three_data);zonghedefen_three = 1:31;  %生成记录点
for i = 1:31%利用（C = D-/((D-)+(D+))）topsis模型与熵权进行结合D_jia  = sqrt(W(1)*(TOP_three_data(i,1)-D_max(1))^2 + W(2)*(TOP_three_data(i,2)-D_max(2))^2 + ...W(3)*(TOP_three_data(i,3)-D_max(3))^2 + W(4)*(TOP_three_data(i,4)-D_max(4))^2 + ...W(5)*(TOP_three_data(i,5)-D_max(5))^2);D_jian = sqrt(W(1)*(TOP_three_data(i,1)-D_min(1))^2 + W(2)*(TOP_three_data(i,2)-D_min(2))^2 + ...W(3)*(TOP_three_data(i,3)-D_min(3))^2 + W(4)*(TOP_three_data(i,4)-D_min(4))^2 + ...W(5)*(TOP_three_data(i,5)-D_min(5))^2);zonghedefen_three(i) = D_jian/(D_jia+D_jian);  % 一次就是产生31个数据
endzonghedefen_three = zonghedefen_three';% 指标四
Z = TOP_four_data;
[n,m] = size(Z);
D = zeros(1,m);  % 初始化保存信息效用值的行向量
for i = 1:mx = Z(:,i);  % 取出第i列的指标p = x / sum(x);% 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以要判断n = length(p);   % 向量的长度lnp = zeros(n,1);   % 初始化最后的结果for j = 1:n   % 开始循环if p(j) == 0   % 如果第i个元素为0lnp(j) = 0;  % 那么返回的第i个结果也为0elselnp(j) = log(p(j));endend%到这里结束e = -sum(p .* lnp) / log(n); % 计算信息熵D(i) = 1- e; % 计算信息效用值
end
W = D ./ sum(D);  % 将信息效用值归一化，得到权重
WD_min = min(TOP_four_data);
D_max = max(TOP_four_data);zonghedefen_four = 1:31;  %生成记录点
for i = 1:31%利用（C = D-/((D-)+(D+))）topsis模型与熵权进行结合D_jia  = sqrt(W(1)*(TOP_four_data(i,1)-D_max(1))^2 + W(2)*(TOP_four_data(i,2)-D_max(2))^2);D_jian = sqrt(W(1)*(TOP_four_data(i,1)-D_min(1))^2 + W(2)*(TOP_four_data(i,2)-D_min(2))^2);zonghedefen_four(i) = D_jian/(D_jia+D_jian);  % 一次就是产生31个数据
endzonghedefen_four = zonghedefen_four'; % 指标五
Z = TOP_five_data;
[n,m] = size(Z);
D = zeros(1,m);  % 初始化保存信息效用值的行向量
for i = 1:mx = Z(:,i);  % 取出第i列的指标p = x / sum(x);% 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以要判断n = length(p);   % 向量的长度lnp = zeros(n,1);   % 初始化最后的结果for j = 1:n   % 开始循环if p(j) == 0   % 如果第i个元素为0lnp(j) = 0;  % 那么返回的第i个结果也为0elselnp(j) = log(p(j));endend%到这里结束e = -sum(p .* lnp) / log(n); % 计算信息熵D(i) = 1- e; % 计算信息效用值
end
W = D ./ sum(D);  % 将信息效用值归一化，得到权重
WD_min = min(TOP_five_data);
D_max = max(TOP_five_data);zonghedefen_five = 1:31;  %生成记录点
for i = 1:31%利用（C = D-/((D-)+(D+))）topsis模型与熵权进行结合D_jia  = sqrt(W(1)*(TOP_five_data(i,1)-D_max(1))^2 + W(2)*(TOP_five_data(i,2)-D_max(2))^2 + ...W(3)*(TOP_five_data(i,3)-D_max(3))^2 + W(4)*(TOP_five_data(i,4)-D_max(4))^2);D_jian = sqrt(W(1)*(TOP_five_data(i,1)-D_min(1))^2 + W(2)*(TOP_five_data(i,2)-D_min(2))^2 + ...W(3)*(TOP_five_data(i,3)-D_min(3))^2 + W(4)*(TOP_five_data(i,4)-D_min(4))^2);zonghedefen_five(i) = D_jian/(D_jia+D_jian);  % 一次就是产生31个数据
endzonghedefen_five = zonghedefen_five'; % 指标六
Z = TOP_six_data;
[n,m] = size(Z);
D = zeros(1,m);  % 初始化保存信息效用值的行向量
for i = 1:mx = Z(:,i);  % 取出第i列的指标p = x / sum(x);% 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以要判断n = length(p);   % 向量的长度lnp = zeros(n,1);   % 初始化最后的结果for j = 1:n   % 开始循环if p(j) == 0   % 如果第i个元素为0lnp(j) = 0;  % 那么返回的第i个结果也为0elselnp(j) = log(p(j));endend%到这里结束e = -sum(p .* lnp) / log(n); % 计算信息熵D(i) = 1- e; % 计算信息效用值
end
W = D ./ sum(D);  % 将信息效用值归一化，得到权重
WD_min = min(TOP_six_data);
D_max = max(TOP_six_data);zonghedefen_six = 1:31;  %生成记录点
for i = 1:31%利用（C = D-/((D-)+(D+))）topsis模型与熵权进行结合D_jia  = sqrt(W(1)*(TOP_six_data(i,1)-D_max(1))^2 + W(2)*(TOP_six_data(i,2)-D_max(2))^2 + ...W(3)*(TOP_six_data(i,3)-D_max(3))^2 + W(4)*(TOP_six_data(i,4)-D_max(4))^2);D_jian = sqrt(W(1)*(TOP_six_data(i,1)-D_min(1))^2 + W(2)*(TOP_six_data(i,2)-D_min(2))^2 + ...W(3)*(TOP_six_data(i,3)-D_min(3))^2 + W(4)*(TOP_six_data(i,4)-D_min(4))^2);zonghedefen_six(i) = D_jian/(D_jia+D_jian);  % 一次就是产生31个数据
endzonghedefen_six = zonghedefen_six';%% 数据整合
zonghedefen = [zonghedefen_one,zonghedefen_two,zonghedefen_three,zonghedefen_four,zonghedefen_five,zonghedefen_six];
% xlswrite('C:\Users\19468\Desktop\writing_two.xlsx',zonghedefen);%% 对一级指标熵权+TOPSIS，保留各一级指标的熵值，信息效用值以及权重，以及最终得分
Z = zonghedefen;
[n,m] = size(Z);
D = zeros(1,m);  % 初始化保存信息效用值的行向量
for i = 1:mx = Z(:,i);  % 取出第i列的指标p = x / sum(x);% 注意，p有可能为0，此时计算ln(p)*p时，Matlab会返回NaN，所以要判断n = length(p);   % 向量的长度lnp = zeros(n,1);   % 初始化最后的结果for j = 1:n   % 开始循环if p(j) == 0   % 如果第i个元素为0lnp(j) = 0;  % 那么返回的第i个结果也为0elselnp(j) = log(p(j));endend%到这里结束e(i) = -sum(p .* lnp) / log(n); % 计算信息熵D(i) = 1- e(i); % 计算信息效用值
end
W = D ./ sum(D);  % 将信息效用值归一化，得到权重
WD_min = min(zonghedefen);
D_max = max(zonghedefen);zonghedefen_last = 1:31;  %生成记录点
for i = 1:31%利用（C = D-/((D-)+(D+))）topsis模型与熵权进行结合D_jia  = sqrt(W(1)*(zonghedefen(i,1)-D_max(1))^2 + W(2)*(zonghedefen(i,2)-D_max(2))^2 + ...W(3)*(zonghedefen(i,3)-D_max(3))^2 + W(4)*(zonghedefen(i,4)-D_max(4))^2 + ...W(5)*(zonghedefen(i,5)-D_max(5))^2 + W(6)*(zonghedefen(i,6)-D_max(6))^2);D_jian = sqrt(W(1)*(zonghedefen(i,1)-D_min(1))^2 + W(2)*(zonghedefen(i,2)-D_min(2))^2 + ...W(3)*(zonghedefen(i,3)-D_min(3))^2 + W(4)*(zonghedefen(i,4)-D_min(4))^2 + ...W(5)*(zonghedefen(i,5)-D_min(5))^2 + W(6)*(zonghedefen(i,6)-D_min(6))^2);zonghedefen_last(i) = D_jian/(D_jia+D_jian);  % 一次就是产生31个数据
endzonghedefen_last = zonghedefen_last'; 
e = e';
D = D';
W = W';%% 肘部法制聚类（一级）
[n, p] = size(zonghedefen_last);
K = 8;
D = zeros(K, 2);
for k = 2:K[label, c, sumd, d] = kmeans(zonghedefen_last, k, 'dist', 'sqeuclidean');% data，n×p原始数据向量% label，n×1向量，聚类结果标签；% c，k×p向量，k个聚类质心的位置% sumd，k×1向量，类间所有点与该类质心点距离之和% d，n×k向量，每个点与聚类质心的距离sse1 = sum(sumd.^2);D(k, 1) = k;D(k, 2) = sse1;
end
% 绘制聚类偏差图
plot(D(2:end, 1), D(2:end, 2), 'o-');
title('不同K值聚类偏差图')
xlabel('分类数(K值)')
ylabel('簇内误差平方和')%% 肘部法制聚类（二级）
[n, p] = size(zonghedefen);
K = 8;
D = zeros(K, 2);
for k = 2:K[label, c, sumd, d] = kmeans(zonghedefen, k, 'dist', 'sqeuclidean');% data，n×p原始数据向量% label，n×1向量，聚类结果标签；% c，k×p向量，k个聚类质心的位置% sumd，k×1向量，类间所有点与该类质心点距离之和% d，n×k向量，每个点与聚类质心的距离sse1 = sum(sumd.^2);D(k, 1) = k;D(k, 2) = sse1;
end
% 绘制聚类偏差图
plot(D(2:end, 1), D(2:end, 2), 'o-');
title('不同K值聚类偏差图')
xlabel('分类数(K值)')
ylabel('簇内误差平方和')%% 系统聚类（一级）
% 执行系统聚类
Z = linkage(zonghedefen_last, 'ward');% 绘制树状图
figure;
dendrogram(Z);
title('Hierarchical Clustering Dendrogram');% 如果你想要裁剪树状图以显示特定数量的聚类，可以使用以下代码
% 例如，我们想要将数据聚类为4个类别
k = 4;
figure;
dendrogram(Z, k);
title(['Hierarchical Clustering Dendrogram with ', num2str(k), ' clusters']);% 裁剪树状图并获取聚类索引
idx = cluster(Z, 'maxclust', k);% 绘制裁剪后的聚类结果
figure;
gscatter(zonghedefen_last(:,1), idx);
title(['Hierarchical Clustering with ', num2str(k), ' clusters']);
xlabel('Feature 1');
ylabel('Feature 2');
legend('Location', 'best');%% 系统聚类（二级）% % 执行系统聚类
% Z = linkage(zonghedefen, 'ward');
% 
% % 绘制树状图
% figure;
% dendrogram(Z);
% title('Hierarchical Clustering Dendrogram');
% 
% % 如果你想要裁剪树状图以显示特定数量的聚类，可以使用以下代码
% % 例如，我们想要将数据聚类为4个类别
% k = 4;
% figure;
% dendrogram(Z, k);
% title(['Hierarchical Clustering Dendrogram with ', num2str(k), ' clusters']);
% 
% % 裁剪树状图并获取聚类索引
% idx = cluster(Z, 'maxclust', k);
% 
% % 绘制裁剪后的聚类结果
% figure;
% gscatter(zonghedefen(:,1),zonghedefen(:,2),zonghedefen(:,3),zonghedefen(:,4),zonghedefen(:,5),zonghedefen(:,6), idx);
% title(['Hierarchical Clustering with ', num2str(k), ' clusters']);
% xlabel('Feature 1');
% ylabel('Feature 2');
% legend('Location', 'best');