文章目录
前言
谈一下自己为什么要参加MathorCup比赛
第一次参加数学建模比赛是深圳杯练手，当时和韦哥，马总一起为了国赛练手，毫无任何建模知识的人选择了比较难的与机器学习相关的“羊犬逃逸”问题。
之后参加建模国赛（我作为建模手），选择的是数据处理方面的优化问题，最终通过线性规划和多目标规划等模型解出了比较靠谱的优化方案，获得了省一。
上学期（大三上），选修了数学建模和MATLAB两门数学建模的基础课程，才更加的了解到自己拿到的一个实际的应用问题产生的灵感怎么用具体的方法和技巧去实现。
我也算是经过“机缘巧合”了解到MathorCup比赛，高校建模比赛（“小国赛”）。作为一个熟悉点建模，熟悉点编程，主攻论文的热爱数学建模的人，还是耐不住“手痒”，凑齐队友，开始比赛。这次我选择的是D题（18万数据的处理题）。
以下为参赛论文的小部分重述，不对之处还请大佬在评论区或者私信我批评指正！
一、摘要
本题主要研究基站选址规划和区域聚类问题。综合考虑减少成本、提高网络覆盖业务量等要求，增设新基站并合理布局满足移动网络通信需求。通过某实际覆盖网络区域的182807个弱覆盖点和1474个现有基站的数据分析与聚类，构建遗传算法模型，制定相应的目标体系和约束条件，调整相关参数，使用 MATLAB、RStudio等软件求解，得到该区域内新站址的最优规划方案。
二、问题
在一个实际网络中，由于通信带宽增大导致基站覆盖范围减小，在给定区域内出现许多弱覆盖区域，要求在业务量优先、成本减少条件下规划建设新基站解决该网络中的弱覆盖问题。具体问题如下：
1、附件1给出确定区域中的弱覆盖点的坐标和业务量信息，附件2给出区域内已有基站的坐标信息，考虑基站类型和站址之间距离门限条件，在给定区域进行站址规划，达到规划基站解决若覆盖点业务量满足总业务量的90%。
2、结合实际基站覆盖状况，使基站由圆形覆盖变成扇形覆盖，扇形覆盖包含三个扇区及对应扇区主方向，在满足扇区间角度约束和扇区覆盖递减条件下，结合问题一条件，给出最优站址、扇形角度、最大解决弱覆盖点的总业务量比例。
3、使用区域聚类技术将所有弱覆盖点进行聚类，划分多个弱覆盖区域，并降低方法的总时间复杂度。
三、聚类遗传模型
1.模型流程分析
2.数据分析
182807个弱覆盖点和1474个现有站点绘制如图。可以看到弱覆盖点的分步呈现一定的聚集性，所以聚类算法是可行的，而且数据点的密集程度高且不均匀。
3.DBSCAN聚类算法
-
参数说明：
（1）minPoints—最小圈住点的个数，单弱覆盖点（距离其他若覆盖点距离很远）我不认为它是一个异常点。
（2）epsilon—聚类半径，即聚类内弱覆盖点间距离阈值。宏基站的覆盖范围为30，微基站覆盖范围为10，聚类之后的每个类中考虑规划新建基站坐标和基站类型，可以在每个类中建多个微基站与宏基站，所以聚类半径设置为宏基站范围，保证聚类不影响基站类型的选择。 -
DBSCAN聚类MATLAB代码
clear;clc;
%导入弱覆盖点数据表
weak = csvread('附件1 弱覆盖栅格数据(筛选).csv', 1, 0);%行偏移量为1，偏移量为0
plot(weak(:, 1), weak(:, 2), 'b.');
hold on
%导入现网站址数据表
exist = csvread('附件2 现网站址坐标(筛选).csv', 1, 0);%行偏移量为1，偏移量为0
plot(exist(:, 2), exist(:, 3), 'g.');
legend('弱覆盖点', '现网站址')
xlabel('x');ylabel('y');
title('数据分析');
saveas(gcf,'问题一_弱点覆盖及现有站点数据图.png');
grid on
axis equal
axis([0 2500 0 2500])
%DBscan聚类实现
%基于密度的噪声应用空间聚类(DBSCAN)是一种无监督的ML聚类算法。
%无监督的意思是它不使用预先标记的目标来聚类数据点。
position=[weak(:, 1), weak(:, 2)];%需要聚类的弱覆盖点坐标数据
flag = zeros(length(weak),1);%聚类访问标记，表示对应弱覆盖点类别
clusterNum = 0;%聚类出的类的个数（序号）
epsilon = 30;%聚类（社区）的最大半径，30—宏基站的覆盖范围
for i=1:length(position)
pointNow = position(i,:);%初始化聚类半径内的邻域点队列
if flag(i)==0
clusterNum = clusterNum+1;%类的个数(也是类的序号)加1
pStart = 1;%队列起始指针
preFlag = flag;%聚类标记更新
while pStart<=length(pointNow) %判断是否完成队列遍历
curp = pointNow(pStart,:); %得到当前要处理的点
pStart = pStart+1; %队列指针更新
diffp = position-curp; %直接和所有数据比较
dis = sqrt(diffp(:,1).*diffp(:,1)+diffp(:,2).*diffp(:,2));%判断当前点与所有点之间的距离
ind = dis<epsilon; %得到距离小于阈值的索引
flag(ind) = clusterNum; %设置当前聚类标记
diff_ind = (preFlag-flag)<0;%判断本次循环相比上次循环增加的点
tmp = zeros(length(position),1);
tmp(diff_ind) = clusterNum;
flag = flag + tmp; %增加点将其标记为一类
preFlag = flag; %聚类标记更新
pointNow = [pointNow;position(diff_ind,:)];%增加聚类半径内的邻域点队列
end
end
end
%运行时间:20分钟左右，用save命令保存聚类结果
save('flag.mat','flag');%双击flag.mat加载变量到工作区
%clusterNum=357
%DBscan聚类结果
figure;
for i = 1: clusterNum
color = [rand(),rand(),rand()];
data_class = position(find(flag==i),:);
plot(data_class(:,1),data_class(:,2),'.','Color',color,'MarkerFaceColor',color);
hold on
end
plot(data_class(:,1),data_class(:,2),'k*');
hold on
grid on
daspect([1 1 1]);
xlabel('x');ylabel('y');
title('DBscan结果');
%保存当前聚类图窗
saveas(gcf,'DBscan30.png');
- 聚类结果
4.遗传算法
1.MATLAB 代码实现
在这里贴出核心代码，均为原创，具体遗传算法的参数部分隐去，望大家理解
%在每个类中做具体算法处理，将处理结果存入规划站点表
for i = 1:clusterNum
count=count+1;
%clear; clc
%i=1;%单次测试
%清除变量big, bigindex, bigpos, bigW, class_w, class_weak, class_e, class_exist, fitvalue, ms, newpop, pop, popx1, popy1, resultCType, resultW, small, smallindex, smallpos, smallW, weak, x_inrange, xypos, xypos_result, xypos_Result, y_inrange
clearvars -except weak exist flag Plan ResultIndex i clusterNum count;%保留变量
class_weak=weak(flag==i,:);%取出一个类的弱覆盖点
class_xmin=min(class_weak(:,1));
class_xmax=max(class_weak(:,1));
class_ymin=min(class_weak(:,2));
class_ymax=max(class_weak(:,2));
x_inrange=find((exist(:,2)>class_xmin)&(exist(:,2)<class_xmax));
y_inrange=find((exist(:,3)>class_ymin)&(exist(:,3)<class_ymax));
xy_inrange=intersect(x_inrange,y_inrange);
class_exist=exist(xy_inrange,2:3);
%plot(class_weak(:, 1), class_weak(:, 2), 'r.');
%hold on
%plot(class_exist(:, 1), class_exist(:, 2), 'g.');
%遗传算法
%当前类中存在的现有基站数据表，平移到坐标原点，减少染色体长度
class_e(:,1)=class_exist(:,1)-class_xmin;
class_e(:,2)=class_exist(:,2)-class_ymin;
%当前类中存在的弱覆盖点数据表，平移到坐标原点，减少染色体长度
class_w(:,1)=class_weak(:,1)-class_xmin;
class_w(:,2)=class_weak(:,2)-class_ymin;
class_w(:,3)=class_weak(:,3);
%plot(class_w(:, 1), class_w(:, 2), 'r.');
%hold on
%plot(class_e(:, 1), class_e(:, 2), 'g.');
%初始化
rangelength=max([class_xmax-class_xmin,class_ymax-class_ymin]);%xy最长的相对距离差
popsize=;%群体大小
G=;%迭代次数
pc=;%交叉概率
pm=;%变异概率(群体大小不够，变异调高保证群体优化)
xylen=length(dec2bin(rangelength));%xy需要的最长的表示x或y二值字符串位数
chromlength=1+2*xylen;%二值数长度（染色体长度）
%class_xmax;class_xmin
%class_ymax;class_xmax
pop=round(rand(popsize,chromlength));%rand随机产生每个单元为{0,1},行数为popsize，列数为chromlength的矩阵
%遗传主程序
for k=1:G
%*************第一步：评估种群中个体适应度*************
%存最后一位0/1（基站类型）
resultCType=pop(:,end);
%区分最后一位0/1
bigindex=find(pop(:,end)==1);%最后一位为1的个体，1表示覆盖范围30
smallindex=find(pop(:,end)==0);%最后一位为0的个体，0表示覆盖范围10
%取每个pop值，前xylen位转十进制x，接着1xylen位转十进制y
for a1=1:xylen
popx1(:,a1)=2.^(xylen-a1).*pop(:,a1);
end
xypos(:,1)=sum(popx1,2);%按行求和得出x
for a2=xylen+1:chromlength-1
popy1(:,a2)=2.^(chromlength-1-a2).*pop(:,a2);
end
xypos(:,2)=sum(popy1,2);%按行求和得出y
%plot(xypos(:,1),xypos(:,2),'b*');
%计算距离小于等于10，30范围业务量加和W
%(1)第一个指标，单位面积的W
bigpos=xypos(bigindex,:);
smallpos=xypos(smallindex,:);
bigW=zeros(length(bigpos(:,1)),1);
pos1=1;
if ~isempty(class_w) && ~isempty(bigpos)
for a1=1:length(bigpos(:,1))
for a2=1:length(class_w(:,1))
if sqrt((class_w(a2,1)-bigpos(a1,1))^2+(class_w(a2,2)-bigpos(a1,2))^2)<30
bigW(pos1)=bigW(pos1)+class_w(a2,3);
end
end
pos1=pos1+1;
end
end
bw_all=bigW;%存每个点总业务量
bigW=bigW/(30*30*pi);
smallW=zeros(length(smallpos(:,1)),1);
pos1=1;
if ~isempty(class_w) && ~isempty(smallpos)
for a1=1:length(smallpos(:,1))
for a2=1:length(class_w(:,1))
if sqrt((class_w(a2,1)-smallpos(a1,1))^2+(class_w(a2,2)-smallpos(a1,2))^2)<10
smallW(pos1)=smallW(pos1)+class_w(a2,3);
end
end
pos1=pos1+1;
end
end
sw_all=smallW;%存每个点总业务量
smallW=smallW/(10*10*pi);
%(2)第二个指标，1/单位面积成本
bigW=bigW*30/10;%修正???
smallW=smallW*10/1;
%(1)+(2)=fitvalue(每一个矩阵)
big=[bigW,bigindex];
small=[smallW,smallindex];
fitvalue=sortrows([big;small],2);
fitvalue=fitvalue(:,1);
fitvalue=fitvalue';%转置成行
bw=[bw_all,bigindex];
sw=[sw_all,smallindex];
resultW=sortrows([bw;sw],2);%存每个点包含的总业务量
%*************第二步：选择复制*************
totalfit=sum(fitvalue);%求适应值之和
fitvalue=fitvalue/totalfit;%单个个体被选择的概率
fitvalue=cumsum(fitvalue);%如fitvalue=[0.1 0.2 0.3 0.4]，则 cumsum(fitvalue)=[0.1 0.3 0.6 1.0]
[px,py]=size(pop);
ms=sort(rand(px,1));%从小到大排列
fitin=1;
newin=1;
while newin<=px
if ms(newin)<fitvalue(fitin)
newpop(newin,:)=pop(fitin,:);
newin=newin+1;
else
fitin=fitin+1;
end
end
%*************第三步：交叉*************
for a=1:2:px-1
if(rand<pc)
cpoint=round(rand*py);
newpop(a,:)=[newpop(a,1:cpoint),newpop(a+1,cpoint+1:py)];
newpop(a+1,:)=[newpop(a+1,1:cpoint),newpop(a,cpoint+1:py)];
else
newpop(a,:)=newpop(a);
newpop(a+1,:)=newpop(a+1);
end
end
%*************第四步：变异*************
for a=1:px
if(rand<pm)
mpoint=round(rand*py);%确定变异点
if mpoint<=0
mpoint=1;
end
newpop(a,mpoint)=~newpop(a,mpoint);
end
end
%*************进行下一次迭代*************
pop=newpop;
end
%存入待选规划基站表
%取每个pop值，前xylen位转十进制x，接着1xylen位转十进制y
for a1=1:xylen
popx1(:,a1)=2.^(xylen-a1).*pop(:,a1);
end
xypos_result(:,1)=sum(popx1,2);%按行求和得出x
for a2=xylen+1:chromlength-1
popy1(:,a2)=2.^(chromlength-1-a2).*pop(:,a2);
end
xypos_result(:,2)=sum(popy1,2);%按行求和得出y
xypos_result(:,1)=xypos_result(:,1)+class_xmin;
xypos_result(:,2)=xypos_result(:,2)+class_ymin;
%共5列，第一列为x坐标，第二列为y坐标，第三列为该点覆盖业务量，第四列为顺序索引，第五列为建站类型
xypos_Result=sortrows([xypos_result,resultW,resultCType]);
% xtmp=xypos_result(1,1);
% ytmp=xypos_result(1,2);
% Plan(ResultIndex,1:4)=xypos_Result(1,[1 2 3 5]);
% ResultIndex=ResultIndex+1;%规划表索引加1
xtmp=0;
ytmp=0;
for j=1:length(resultCType)
if (xypos_result(j,1)==xtmp)&&(xypos_result(j,2)==ytmp)
%重复了不添加
else
tmpR=xypos_Result(j,[1 2 3 5]);
if isempty(find(Plan(:,1)==tmpR(1))) && ~(tmpR(3)==0) && (tmpR(1)<=2500) && tmpR(2)<=2500 %业务量为0的点不入
Plan(ResultIndex,:)=tmpR;
xtmp=tmpR(1);%更新标记
ytmp=tmpR(2);%更新标记
ResultIndex=ResultIndex+1;%规划表索引加1
end
tmpR=xypos_Result(j,[1 2 3 5]);
if isempty(find(Plan(:,1)==tmpR(1))) && isempty(class_exist) && ~(tmpR(3)==0) && (tmpR(1)<=2500) && tmpR(2)<=2500 %业务量为0的点不入,且要在范围内
Plan(ResultIndex,:)=tmpR;
xtmp=tmpR(1);%更新标记
ytmp=tmpR(2);%更新标记
ResultIndex=ResultIndex+1;%规划表索引加1
else
ff=0;
if ~isempty(class_exist)%考虑门限
for k=1:length(class_exist(:,1))
if (sqrt((class_exist(k,1)-tmpR(1))^2+(class_exist(k,2)-tmpR(2))^2))>10
ff=1;
end
end
if (ff==0) && ~(tmpR(3)==0) && (tmpR(1)<=2500) && tmpR(2)<=2500 %业务量为0的点不入
Plan(ResultIndex,:)=tmpR;
xtmp=tmpR(1);%更新标记
ytmp=tmpR(2);%更新标记
ResultIndex=ResultIndex+1;%规划表索引加1
end
end
end
end
end
end
- 遗传结果
总结
下面从几个方面来总结一下自己比赛结束的心得。
1.论文能力
从第一次参加数学建模比赛开始，每次的论文基本最终都是由我去主要整合和陈述，平时的大大小小的课程设计小组，我也经常在报告和论文中承担比较重要的角色。
我知道很多同学或者说和我有相同学习情景的朋友，不喜欢或者不擅长论文或报告的撰写，但是他们的代码水平要比我好很多。我并不是说自己论文水平有多么好，毕竟咱也还不是研究生，那比起每天搞论文的大佬还是天壤之别，但我对写报告论文至少有比较浓厚的兴趣，和大家分享一下我自己的经验。
论文能力，文书能力这样的写作能力，往往在很多事情的书面沟通和细节处理上起关键的作用。我自己是在学生会的三年里有深切感受倒文书工作的重要性，很多时候要发布的新闻稿或者文案，都是要给全院，全校甚至更多人看的，并不是只有自己明白自己的意思就可以，这需要积累，也需要细心。除了工作，在计算机这个行业，有必要，也必须有一定的论文能力，你需要把自己的编程思想陈述明白，把自己的代码写成规范性的文档。才会有更多人明白你做了什么。
2.抗压能力
其实这次比赛给我感触最深的是，我发现自己的抗压能力显著提高了。抗压能力，对于每一个生活着的人，极其重要。我这里仅仅谈作为计算机专业学生在抗压的经验和处理。
很多事情都是残酷的，天赋，就是残酷的。学计算机的人可能在大学感触就会很深，有些人写代码，写算法，就是一遍过，没有bug；但有些人就是简单的代码，并且这些人也是了解算法的思想，但就是写不对代码，bug很多（用一个bug去解决另一个bug）。当然这个是可以通过大量的刷题去弥补的，但是要明白一个道理。
抗压能力在这个时候就很重要，如果你能承受住自己内心对自己的不断怀疑、冲击，那么你的抗压能力足够优秀。
这次比赛我感受到自己做不出来D题，因为确实题目数据量太大，而且我的队友对于编程的熟练度不高，我需要一个人去对抗核心的算法，并且承受自己内心的怀疑，但我从第一天晚上、到第二天晚上、再到第三天晚上（比赛一共四天），我越来越坚定，也越来越放松。
这是我自己总结的，与大家共勉。
3.身体健康
身体健康很重要！！！这次比赛我就拖着不舒服的身体去持续了四天的高强度比赛，产出的效率肯定比现在我身体健康状况下写博客的速度要慢许多。
对于我们可能未来十年二十年都会面临着高强度的科研、敲代码等高强度的工作，保护好自己身体健康是很重要的，没有身体健康，你没有了去写出好的论文，搞出好的算法的本钱。
结尾感谢：建模队友（爱吃桃子、蒙眼的兔子）
