【比赛总结】从编程位队长的角度看待第十三届华中杯数学建模比赛A题

前言

有幸以编程位和队长的身份大一就参加了一次数学建模比赛，这次比赛是"华中杯"，所以第一次打还是比较有新鲜感和有很多收获的，故记于此。

因为……在前期找指导老师的时候一说是大一的队伍就不建议我们参加（无语ing

所以我们队（至少我）在怀疑我们这么努力会不会有成绩，但是事实是，我们都成长和学到了好多

代码会放在附录位置

本文还会以队长的身份对以后建模提一些自己的看法和意见

我们这次选择的A题（具体题目可以搜到，今年开始华中杯开始公布题目）

（最后回来写：已经写完了，就先写这些吧）

我们的做法

A题是一个关于颜色RGB值的问题，刚拿到问题没仔细分析的时候我误判为图论

（虽然后来看群里有人讨论拓扑结构，还听到别的队有讨论最短路的，可以说图论不一定是误判哦）

对于A题的第一问，我们开始考虑直接用RGB空间坐标系中的欧氏距离表示颜色的相似度

后经过对比误差略大，后来我们采取了推导了一个线性表达式并给相关系数求出来的结果使误差减小

编程位可以干什么呢？可以看到第一问的关键在于如何定义颜色的相似度，编程位的作用就是

把你们队想到的计算相似度的公式用代码语言跑出来，其实这一点稍微有点语言基础基本都可以做到

有一些我们没有实现的算法，我们没有想到但是可以去猜测一种做法：某个队根据几个值得出了一个非线性公式

并有多个限制条件（非线性规划）然后配合模拟退火等智能优化算法得出更优解

只可惜我们队的编程位水平太低（没错就是我）不能很好的找到公式和完美运用模拟退火……

对于A题的第二问，我们还讨论了好久这一问究竟是10小问还是11小问……

我觉得第二题是整个A题的重点，应该是编程和算法的关键，也是评分的关键

说一下编程的角度看待这个问题：

可以很容易想到，把22种颜色分类，去找1600万个点中距离这些类的“重心”最远的那个点一定是第一个增加的

这里所谓的“重心”就是建立合适的坐标系求得的中心

我们队建立的坐标系是：RGB坐标系->HSV坐标系->圆锥坐标系下HSV坐标系（有理由的）

其实关键在于如何分类。开始想到了神经网络，奈何我的水平太差写不出来

后来发现kmeans算法我可以很快的实现，就用kmeans去分类

分好类之后我们可以对于每一个类去二分球的半径，不在所有球内的就是答案

至于如何判断是否在求内，其实O(n)跑一遍所有颜色到22种颜色的重心排序后的结果就好

下面是最后论文里我提供给论文位的step，可以看看吧。

step1.把22种瓷砖颜色的RGB值转化为圆锥坐标系下坐标step2.枚举类数K∈[1,22]，用Kmeans算法求得重心坐标step3.计算22种已有颜色到重心的距离和，并确定合适的K值step4.枚举增加一个其他颜色的RGB值，用Kmeans算法求得重心坐标step5.计算重心偏移量，降序排序并把偏移量最大颜色设为已有（必选）step6.二分球的半径并判断其他颜色是否在已有颜色画出的球内step7.若不在则设为已有（必选），重复操作6，直到球的半径在误差精度内

对于A题的第三问，我们在第二问的基础上给出了11种颜色，并（基本可以）确定第三问的答案一定在这11种中。（至于为什么不是第二问题目问的10种，也是有原因的哦）

然后用二进制去表示选还是不选，一共有2047种状态（即2^11-1，全是0的状态也不可)，美其名曰0-1规划罢了

然后把这些状态跑一遍，用权重和除以1的个数即为该状态下综合权值（综合考虑成本和表达效果）

具体解释上面这句话：

权重和：权重哪里来的？我们用了比较主观的层次分析法（为什么用主观点的也需要解释理由哦）

用层次分析法分析H、S、V三者在决定颜色中的权重，并对于每一个因素又分析了11种颜色的权重，并分别进行一致性检验

最后给出每一种颜色的权重

值得一提的是，第一问我们直接用的RGB，后面两问我们采取了HSV（为什么这样选其实就蕴含在二者的对比之中）

而权重和就是每个状态二进制是1的颜色权重相加

因为每种颜色的成本一样，那么总状态的成本可以用颜色的个数，即1的个数来表示

这里查到一个类似的定义“投资回报率”，“收入/成本”可以相当于我们的“权重和/1的个数”

扯了这么多思路，那么编程位干什么呢？O(nlogn)排序（其实O(n^2)也没事，毕竟三天呢肯定能跑出来）

这里好像学过语言的好像也都能写了哈哈哈

然后取出最后总权重最大的那个状态对应的颜色就好了

学到/成长了什么

这一部分会再开一篇博文去详细解释的

第一次领导数模的队伍，对于做队长有一点自己的想法了（下一部分说）

认识了两个有趣且有能力的灵魂

在准备的过程中学到Matlab的好多操作，也学习了遗传算法、粒子群、模拟退火等智能优化算法

比赛的时候深入了解到RGB、HSV等颜色模型

学习了并实现了Kmeans算法，了解了神经网络的原理（虽然现在还不能实现

第一次为了一场比赛睡眠不足，体会到了付出的感觉

跟着论文位学会了查知网、SCI等，并学到了一些写论文的操作、步骤和技巧

学会了一点点excel的操作，会用ProcessOn去画流程图，懂得了markdown的强大哈哈哈

会用markdown画mermaid了

对以后建模队伍的想法与建议

建议编程位不要当队长，建议论文位当

论文位应该善于去寻求其他两个队员的帮助，时不时地提出自己的想法和需要

应该三个人都明白每个题的做法，不明白也要讲明白不然最后的评价推广摘要等不好写

如果有了分歧不要急着否定一定要听取队友的看法，尊重队友的做法和想法

这场比赛编程的源程序

层次分析法 AHP.m

clear; clc;
fid=fopen('mat1.txt','r');
%n1行 n2列 mat1:n1=n2=3 , mat2/mat3/mat4:n1=n2=11
n1=3;n2=3;
A=[];
for i=1:n1 tmp=str2num(fgetl(fid)); A=[A;tmp]; %读准则层判断矩阵
end
[n,n]=size(A);
x=ones(n,100);
y=ones(n,100);
m=zeros(1,100);
m(1)=max(x(:,1));
y(:,1)=x(:,1);
x(:,2)=A*y(:,1);
m(2)=max(x(:,2));
y(:,2)=x(:,2)/m(2);
p=0.0001; i=2; k=abs(m(2)-m(1));
while k>pi=i+1;x(:,i)=A*y(:,i-1);m(i)=max(x(:,i));y(:,i)=x(:,i)/m(i);k=abs(m(i)-m(i-1));
end
a=sum(y(:,i));
w=y(:,i)/a;
t=m(i);
disp(w);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];
CR=CI/RI(n);
if CR<0.10disp('Accept');disp('CI=');disp(CI);disp('CR=');disp(CR);
end

问题1求解程序 problem1.cpp

#include<iostream>
#include<cstdio>
#include<fstream>
#include<cstring>
#include<string>
#include<cmath>
#define INF 2147483647
using namespace std;struct Node {int r , g , b;Node(){};Node(int rr , int gg , int bb) {r=rr; g=gg; b=bb;}
};
Node t[22];
const double d=64;
const double mov = 0.01;
const double k = 0.9;
int n = 0;
string s;void init() {t[0] = Node(0,0,0);t[1] = Node(255,255,255);t[2] = Node(255,0,0);t[3] = Node(246,232,9);t[4] = Node(72,176,64);t[5] = Node(27,115,186);t[6] = Node(53,118,84);t[7] = Node(244,181,208);t[8] = Node(255,145,0);t[9] = Node(177,125,85);t[10] = Node(92,59,144);t[11] = Node(11,222,222);t[12] = Node(228,0,130);t[13] = Node(255,218,32);t[14] = Node(118,238,0);t[15] = Node(17,168,226);t[16] = Node(255,110,0);t[17] = Node(201,202,202);t[18] = Node(255,249,177);t[19] = Node(179,226,242);t[20] = Node(249,225,214);t[21] = Node(186,149,195);
}Node sol(string ss) {int zx = ss.find('(') , len = ss.size()-1 , cnt = 0;int flag1 = 0 , flag2 = 0 , flag3 = 0;Node tmp;for(int i=zx+1; i<=len; i++) {if(ss[i]==','||ss[i]==')') {if(!flag1&&!flag2&&!flag3) {tmp.r = cnt;cnt = 0;flag1 = 1;continue;}if(flag1&&!flag2&&!flag3) {tmp.g = cnt;cnt = 0;flag2 = 1;continue;}if(flag1&&flag2&&!flag3) {tmp.b = cnt;cnt = 0;flag3 = 1;continue;}} else {cnt = (cnt*10)+(ss[i]-'0');}}return tmp;
}int main() {//输入输出数据 ifstream in("附件2：图像1颜色列表.txt",ios::in);ofstream out("result1.txt",ios::out);init(); in>>s;out<<"序号,瓷砖颜色编号"<<endl;while(in>>s) {n++;Node tmp = sol(s);double ans1 = INF , ans2 = INF;int opt;for(int j=0; j<=21; j++) {int dism = max(abs(t[j].r-tmp.r) , max(abs(t[j].g-tmp.g) , abs(t[j].b-tmp.b)));double r1 = tmp.r+mov , r2 = t[j].r+mov;double g1 = tmp.g+mov , g2 = t[j].g+mov;double b1 = tmp.b+mov , b2 = t[j].b+mov;double l1=sqrt(r1*r1+g1*g1+b1*b1);   double l2=sqrt(r2*r2+g2*g2+b2*b2);double disn = sqrt((t[j].r-tmp.r)*(t[j].r-tmp.r)+(t[j].g-tmp.g)*(t[j].g-tmp.g)+(t[j].b-tmp.b)*(t[j].b-tmp.b));double sita = 1.0*(r1*r2+g1*g2+b1*b2)/(l1*l2);if(dism<=d) {if(disn*k+sita*(k-1)<ans1) {ans1 = disn*k+sita*(k-1);ans2 = dism;opt = j+1;} } else {if(dism<ans2) {opt = j+1;ans2 = dism;}}}out<<n<<","<<opt<<endl;  }in.close();out.close();return 0;
}

问题2求解程序 problem2.cpp

#include<iostream>
#include<cstdio>
#include<string>
#include<cmath>
#include<cstdlib>
#include<ctime>
#include<fstream>
#include<map>
#include<queue>
#include<algorithm>
#include"Kmeans.h"
using namespace std;#define N 25
#define K 15
int Kk;
extern Point point[N] , mean[15];
extern int center[N];
int tx , ty , tz; //x:H y:S z:V
int m , vis[260][260][260];
Point ans[N+200];
const double eps = 0.2;struct Node {Point x;double y;Node() {};Node(Point xx , double yy) {x=xx , y=yy;}
};bool operator < (Node a , Node b) {return a.y>b.y;
}
int cp;
Node q[17000000];
map<Point,Point>mp;void Kmeans() {srand((unsigned int)time(NULL));for(int i=0; i<Kk; i++) {int j = rand() % Kk;mean[i].x = point[j].x;mean[i].y = point[j].y;mean[i].z = point[j].z;}double temp1 , temp2;Cluster();temp1 = GetE(); GetMean(center);Cluster();temp2 = GetE();while(fabs(temp2-temp1)!= 0) {temp1 = temp2;GetMean(center);Cluster();temp2 = GetE();}
}int jud(double R , int n , Point a[]) {int cnt = m;a[cnt] = q[0].x;cnt++;for(int i=1; i<cp; i++) {if(q[i].y<=eps) break; int flag = 0;for(int j=0; j<cnt; j++) {if(GetDistance(a[j] , q[i].x)<=R) {flag = 1;break;}}if(!flag) {a[cnt] = q[i].x;cnt++;     }}if(cnt-m>=n)return 1;return 0;
}int main() {ifstream in("data.txt",ios::in); //读入22种颜色的RGB值 ofstream out("out.txt",ios::out); //将该做法的答案输出到out.txt中 while(in>>tx>>ty>>tz) {point[m] = RGBtoHSV(tx,ty,tz);vis[tx][ty][tz] = 1;m++;}Kk=8;Kmeans();Point tep = mean[0];for(int R=0; R<=255; R++)for(int G=0; G<=255; G++)for(int B=0; B<=255; B++) {if(vis[R][G][B]) continue;point[m]=RGBtoHSV(R,G,B);m++;Kmeans();m--;q[cp++] = Node(point[m],sqrt((tep.x-mean[0].x)*(tep.x-mean[0].x)+(tep.y-mean[0].y)*(tep.y-mean[0].y)+(tep.z-mean[0].z)*(tep.z-mean[0].z)));mp[point[m]] = Point(R,G,B);}sort(q,q+cp);for(int i=0; i<m; i++)ans[i] = point[i];for(int i=1; i<=10; i++) {double l=0 , r =255;out<<i<<":"<<endl;while(fabs(r-l)>eps) {double mid = (l+r)/2;if(jud(mid,i,ans))l = mid;elser = mid;}for(int j=m; j<m+i; j++)out<<mp[ans[j]].x<<" "<<mp[ans[j]].y<<" "<<mp[ans[j]].z<<endl;out<<endl;}in.close();out.close();return 0;
}

问题3求解程序 problem3.cpp

#include<iostream>
#include<fstream>
#include<algorithm>
#include<cstring>
using namespace std;int init[15];
double qz1[15] , qz2[15] , qz3[15];
struct Node {int no;double f;
};
Node m[2050];
const double k1 = 0.6483; //H S V各自的权值
const double k2 = 0.1220;
const double k3 = 0.2297;bool operator < (Node& a, Node& b) {return a.f < b.f;
}void sol(int num) {memset(init , 0 , sizeof(init)); int cnt = 0;while (num>=2) {init[cnt++] = num % 2;num >>= 1;}init[cnt++] = num;
}int main() {ifstream in("weight.txt", ios::in);for (int i = 0; i < 11; i++)in>>qz1[i]>>qz2[i]>>qz3[i];for (int i = 1; i <= 2047; i++) {sol(i);int n = 0;double sum = 0;for (int j = 0; j < 11; j++) {if (init[j] == 1) {n++;sum += qz1[j]*k1+qz2[j]*k2+qz3[j]*k3;}}m[i] = {i , sum/n};}sort(m + 1, m + 2048);for(int i=2047; i>=2044; i--) { //输出L/n最大的前几个 sol(m[i].no);for(int i=0; i<11; i++)cout<<init[i]<<" ";cout<<endl;cout<<m[i].f<<endl<<endl;}in.close();return 0;
}

Kmeans算法头文件 Kmeans.h

#include<iostream>
#include<cstdio>
#include<string>
#include<cmath>
#include<cstdlib>
#include<ctime>
#include<fstream>
using namespace std;#define PI acos(-1)
#define N 25
extern int Kk;
#define K 15struct Point {double x , y , z;Point() {};Point(double xx , double yy , double zz) {x=xx , y=yy , z=zz;}
};bool operator < (Point a , Point b) {return a.x<b.x;
}Point point[N] , mean[15];
int center[N];
extern int m;double GetDistance(Point point1, Point point2) {return sqrt((point1.x-point2.x)*(point1.x-point2.x)+(point1.y-point2.y)*(point1.y-point2.y)+(point1.z-point2.z)*(point1.z-point2.z));
}void GetMean(int center[N]) { // 计算每个簇的中心点Point tmp;for(int i=0; i<Kk; i++) {int count = 0; tmp.x = tmp.y = tmp.z = 0.0;for(int j=0; j<m; j++) {if(i==center[j]) {count++;tmp.x += point[j].x;tmp.y += point[j].y;tmp.z += point[j].z;}}tmp.x /= count; tmp.y /= count; tmp.z /= count;mean[i] = tmp;}
}double GetE() { // 计算平方误差函数double cnt = 0.0, sum = 0.0;for(int i=0; i<Kk; i++) {for(int j=0; j<m; j++) {if(i==center[j]) {cnt=(point[j].x-mean[i].x)*(point[j].x-mean[i].x)+(point[j].y-mean[i].y)*(point[j].y-mean[i].y)+(point[j].z-mean[i].z)*(point[j].z-mean[i].z);sum += cnt;}}}return sum;
}void Cluster() { // 把N个点聚类double minn;double distance[N][K];for(int i=0; i<m; i++) {minn = 999999.0;for(int j=0; j<Kk; j++)distance[i][j] = GetDistance(point[i], mean[j]);for(int j=0; j<Kk; j++) {if(distance[i][j] < minn) {minn = distance[i][j];center[i] = j;}}}
}const double RR = 100;
const double angle = 30;
double h = RR*cos(angle/180*PI) , r = RR*sin(angle/180*PI); //圆锥坐标系的参数 double mymax(double a , double b , double c) {return max(a , max(b , c));
}double mymin(double a , double b , double c) {return min(a , min(b , c));
}Point RGBtoHSV(double R , double G , double B) { //RGB转HSV圆锥坐标系下坐标 R/=255; G/=255; B/=255;double H , S , V;double maxx=mymax(R,G,B) , minn=mymin(R,G,B);if(maxx==minn)H = 0;else {if(R==maxx)H = (G-B)/(maxx-minn);if(G==maxx)H = 2+(B-R)/(maxx-minn);if(B==maxx)H = 4+(R-G)/(maxx-minn);H = H*60;if(H<0) H+=360; }if(maxx==0)S = 0;elseS = (maxx-minn)/maxx;V = maxx;return Point(r*V*S*cos(H),r*V*S*sin(H),h*(1-V));
}