2019年全美数学建模竞赛（B题，重心法选址+类蚁群算法路径规划）

整体论文研究了半天没法上传，就把中间的几部分的代码上传了。

1.无人机属性预处理

#include<bits/stdc++.h>
using namespace std;
#define maxn 62
struct point{
    int x,y,z,w;
};
struct plane{
    int op;//装载类型
    int x,y,z;//长宽高
    double w,v,t;//载重，速度，飞行时间
    plane(){};
    plane(int tx,int ty,int tz,double tw,double tv,double tt,int top){
        x=tx;y=ty;z=tz;w=tw;v=tv;t=tt;op=top;
    };
    out(){
        printf("(%d,%d,%d)\n",x,y,z);
    }
}a[10];
struct Medcine_pok{
    int x,y,z,w;
    Medcine_pok(){};
    Medcine_pok(int xx,int yy,int zz,int ww){x=xx;y=yy;z=zz;w=ww;}
}mdk[3];
void init(){
    mdk[0]=Medcine_pok(14,7,5,2);
    mdk[1]=Medcine_pok(5,8,5,2);
    mdk[2]=Medcine_pok(12,7,4,3);
    a[1]=plane(45,45,25,3.5,40,35,1);
    a[2]=plane(30,30,22,8,79,40,1);
    a[3]=plane(60,50,30,14,64,35,2);
    a[4]=plane(25,20,25,11,60,18,1);
    a[5]=plane(25,20,27,15,60,15,2);
    a[6]=plane(40,40,25,22,79,24,2);
    a[7]=plane(32,32,17,20,64,16,2);
}
int sx,sy,sz;
bool vis[maxn][maxn][maxn];
int b[maxn][maxn];//当前顶位置
double sw;


bool cla(){//枚举把物品p放入箱子x的方式,查看是否可以放开
    //------------------------初始化------------------------------
    memset(vis,0,sizeof(vis));
    memset(b,0,sizeof(b));
    //memset(pvis,0,sizeof(pvis));
    //------------------------------------------------------------
}
int main(){
    init();
    double ans=0;
    for(int i=0;i<8;i++){
            ans+=a[i].v*a[i].t/60/2;
        printf("%d %.2lf\n",a[i].op,a[i].v*a[i].t/60/2);
    }
    printf("%lf\n",ans/7.0);
    return 0;
}

/* 程序目标： 1.对于空间(x,y,z,w)，将若干种物品放入其中，寻找一种最优的放置方案/判断是否可以放置成功。 2.对于空间(x,y,z,w)，在放置某些确定的物品的前提下，如何放置尽量多的添加其他物品？ 3.对于三倍空间(x,y,z,w)，如何将已知物品分配放置到其中，使得占位尽量小/是否可以放置成功。 */

2.无人机选择代码：

#include<bits/stdc++.h>
using namespace std;
#define maxn 202
#define inf 1000222
#define p2(x) ((x)*(x))
const double R=6378.137;
int n,m;
struct point{
    double x,y;
    point(){};
    point(double tx,double ty){
        x=tx,y=ty;
    }
    void out(){
        printf("(%lf,%lf)\n",x,y);
    }
}p[maxn],hs[10];
struct edge{
    int s,t;
}e[maxn];
double operator -(point a,point b){//计算ab两点距离
    double X1=a.y*3.1415926/180;
    double X2=b.y*3.1415926/180;
    double Y1=a.x*3.1415926/180;
    double Y2=b.x*3.1415926/180;
    return R*acos( cos(Y1)*cos(Y2)*cos(X1-X2)+sin(Y1)*sin(Y2) );
}
point operator +(point a,point b){
    return point(a.x+b.x,a.y+b.y);
}
point operator *(int a,point b){
    return point(a*b.x,a*b.y);
}
double f(point x,point aim){
    return pow(u,x-aim);
}
void input(){
    scanf("%lf%lf",&hs[1].x,&hs[1].y);
    scanf("%lf%lf",&hs[2].x,&hs[2].y);
    scanf("%lf%lf",&hs[3].x,&hs[3].y);
    scanf("%d",&n);
    for(int i=0;i<n;i++){
        scanf("%lf%lf",&p[i].x,&p[i].y);
    }
}
void input_edge(){
    scanf("%d",&m);
    for(int i=0;i<m;i++){
        scanf("%d%d",&e[i].s,&e[i].t);
        //printf("(%d,%d)=%lf\n",e[i].s,e[i].t,p[e[i].s]-p[e[i].t]);
    }
}
int ans[5][2004];//0.01km~20km
bool f(edge x,point m,double r){
    return ((p[x.s]-m)-r)*((p[x.t]-m)-r)<0;
}
//17.972856 -67.159519 18.485500 -65.636833
double fs(point a,point b){
    if(a-b>40)return 0;
    return pow(0.6,(a-b) -2 );
}
int main(){
    freopen("dian.txt","r",stdin);
    //freopen("range_s.txt","w",stdout);
    input();
    input_edge();
    double ans=0;
/*-------------------------------------------------------------*/
    for(int i=0;i<m;i++){
        ans+=p[e[i].s]-p[e[i].t];
    }
    printf("%lf\n",ans);
/*-------------------------------------------------------------*/
    double u=(18.485500-17.972856)/600;
    double v=-(-67.159519+65.636833)/1000;
    for(double i=17.972856- 5*u;i<=18.485500+5*u;i+=u){
        for(double j=-67.159519- 5*v;j<=-65.636833+5*v;j+=v){
            double ans=0;
            for(int k=0;k<n;k++){
                ans+=fs(point(i,j),p[k]);
            }
            printf("%lf ",ans);
        }
        printf("\n");
    }
/*-------------------------------------------------------------*/
    double max1=-inf,max2=-inf,min1=inf,min2=inf;
    for(int i=0;i<n;i++){
        max1=max(max1,p[i].x);
        max2=max(max2,p[i].y);
        min1=min(min1,p[i].x);
        min2=min(min2,p[i].y);
    }
    printf("%lf %lf %lf %lf\n",min1,min2,max1,max2);
/*-------------------------------------------------------------*/
    for(int aim=1;aim<=1;aim++){
        for(double x=-20;x<=20;x+=0.3){
            //printf("%lf\n",x);
            for(double y=-20;y<=20;y+=0.3){
                int sum=0;
                double r=sqrt(x*x+y*y);
                for(int i=0;i<m;i++){
                    if(f(e[i],hs[aim],r))sum++;
                }
                printf("%d ",sum);
            }
            printf("\n");
        }
        for(double r=0;r<100;r+=0.01){
            int sum=0;
            for(int i=0;i<m;i++){
                if(f(e[i],hs[aim],r))sum++;
            }
            printf("%d\n",sum);
        }
    }
    return 0;
}
/* 18.297473 -65.932790 18.335043 -66.093564 18.354452 -66.544844 (18.297473,-65.932790) (18.290914,-66.069934) (18.335043,-66.093564) (18.330872,-66.128470) (18.354452,-66.544844) */

3.重心法选择集装箱安放地点

#include<bits/stdc++.h>
using namespace std;
#define maxn 202
#define p2(x) ((x)*(x))
const double R=6378.137;
double u;
int n;
struct point{
    double x,y;
    point(){};
    point(double tx,double ty){
        x=tx,y=ty;
    }
    void out(){
        printf("(%lf,%lf)\n",x,y);
    }
}p[maxn],hop[10];
double operator -(point a,point b){//计算ab两点距离
    double X1=a.y*3.1415926/180;
    double X2=b.y*3.1415926/180;
    double Y1=a.x*3.1415926/180;
    double Y2=b.x*3.1415926/180;
    return R*acos( cos(Y1)*cos(Y2)*cos(X1-X2)+sin(Y1)*sin(Y2) );
}
point operator +(point a,point b){
    return point(a.x+b.x,a.y+b.y);
}
point operator *(int a,point b){
    return point(a*b.x,a*b.y);
}
double f(point x,point aim){
    return pow(u,x-aim);
}
void input(){
    scanf("%d",&n);
    for(int i=1;i<=5;i++){
        scanf("%lf%lf",&hop[i].x,&hop[i].y);
    }
    for(int i=1;i<=n;i++){
        scanf("%lf%lf",&p[i].x,&p[i].y);
    }
}
int main(){
    freopen("dian.txt","r",stdin);
    freopen("fitness_search.txt","w",stdout);
    input();
    for(u=0.01;u<=1;u+=0.01){
        printf("%.2lf\t",u);
        for(int aim=1;aim<=5;aim++){
            point ans;
            double sum=0;
            ans.x=ans.y=0;
            for(int i=1;i<=n;i++){
                sum+=f(p[i],hop[aim]);
            }
            for(int i=1;i<=5;i++){
                sum+=f(hop[i],hop[aim]);
            }
            for(int i=1;i<=n;i++){
                ans.x+=f(p[i],hop[aim])*p[i].x;
                ans.y+=f(p[i],hop[aim])*p[i].y;
            }
            for(int i=1;i<=5;i++){
                ans.x+=f(hop[i],hop[aim])*hop[i].x;
                ans.y+=f(hop[i],hop[aim])*hop[i].y;
            }
            ans.x/=sum;
            ans.y/=sum;
            printf("(%lf,%lf)\t",ans.x,ans.y);
            printf("%lf\t",hop[aim]-ans);
        }
        printf("\n");
    }
    return 0;
}
/* (18.325725,-65.652346) (18.224337,-66.032753) (18.432155,-66.070619) (18.400684,-66.155618) (18.467167,-66.726353) (18.297473,-65.932790) (18.290914,-66.069934) (18.335043,-66.093564) (18.330872,-66.128470) (18.354452,-66.544844) */

4.类蚁群算法解决路径规划问题

#include<bits/stdc++.h>
using namespace std;
#define maxn 305
#define p2(x) ((x)*(x))
const double inf=1e14;
const double R=6378.137;
const double P_jump=0.05;
const double u1=0.5;
const double u2=0.5;
const double eps=7;
const double eg=0.5;
double sum_len;//总边长度
double max_cost;//花费总时间
vector<int>u[maxn];//每个点的出边
int num[maxn][maxn];//每条边走过的次数
int id[maxn][maxn];//每对点对应的边的编号，-1表示没有
double dis[maxn][maxn];//两点之间的最近距离
double con[maxn][maxn];//两点之间的直接距离
struct edge{
    int s,t;//s-->t
    double w;
    edge(){};
    edge(int ss,int tt,double ww){s=ss;t=tt;w=ww;}
}e[maxn];
int nume;//边数[0,nume)
struct point{
    int id;//该点的id
    double lx,ly,x,y;///纬度x,经度y,图片横坐标lx，纵坐标ly
    point(){};
    point(double tx,double ty){
        x=tx,y=ty;
    }
    void out(){
        printf("(%lf,%lf)\n",x,y);
    }
}p[maxn];
int nump;//点数[0,nump)
struct plane{
    double V,T;//速度,飞行时间
    bool flag;//是否有载重
    int s,t;//起点，目标终点
    double len,res;//能走最长距离len，剩余可以距离res
    double wmax,w;//最大载重量wmax,实际任务载重量w
    vector<edge>way;//已走路径
    point now;//当前所在位置
    double cost(){//当前所用时间
        double ret=0;
        for(int i=0;i<way.size();i++){
            ret+=con[way[i].s][way[i].t];
        }
        return ret;
    }
}ant[maxn];
int numa;//飞机个数
bool operator != (point a,point b){
    return a.x!=b.x||a.y!=b.y;
}
bool operator == (point a,point b){
    return a.x==b.x&&a.y==b.y;
}
double operator -(point a,point b){//计算ab两点距离
    double ret=0;
    double X1=a.y*3.1415926/180;
    double X2=b.y*3.1415926/180;
    double Y1=a.x*3.1415926/180;
    double Y2=b.x*3.1415926/180;
    //ret=sqrt(p2(sin((a.x-b.x)/2.0))+ cos(a.x)*cos(b.x)*p2(sin((a.y-b.y)/2.0)) );
    //ret=2*R*asin(ret);
    return R*acos( cos(Y1)*cos(Y2)*cos(X1-X2)+sin(Y1)*sin(Y2) );
}
void add_edge(int s,int t,double w){
    id[s][t]=id[t][s]=nume;
    con[s][t]=con[t][s]=w;
    e[nume++]=edge(s,t,w);
    u[s].push_back(t);
    u[t].push_back(s);
    sum_len+=w;
}
void input(){
    int T;
    freopen("in2.txt","r",stdin);
    freopen("in_ans.txt","w",stdout);
    scanf("%d%d%d",&nump,&T,&numa);
    for(int i=0;i<nump;i++){//0~4为医院
        scanf("%lf%lf",&p[i].x,&p[i].y);
        p[i].id=i;
    }
    while(T--){//道路
        int ts,tt;  double tw;
        scanf("%d%d",&ts,&tt);
        tw=p[ts]-p[tt];
        printf("%d %d %lf\n",ts,tt,tw);
        add_edge(ts,tt,tw);
    }
    for(int i=0;i<numa;i++){//若干架飞机
        plane temp;
        scanf("%d%d%lf%lf%lf%lf",&temp.s,&temp.t,&temp.V,&temp.T,&temp.wmax,&temp.w);
        temp.flag=1;
        temp.len=temp.V*temp.T/60;
        temp.res=temp.len;
        temp.now=p[temp.s];
        ant[i]=temp;
    }

}
void init(){
    max_cost=0;
    sum_len=0;
    nume=0;//边数为0
    memset(id,-1,sizeof(id));//清空边
    memset(num,0,sizeof(num));//清空走过路径
    for(int i=0;i<maxn;i++){//清空邻接表边
        u[i].clear();
    }
}
void floyd(){
    int n=nump;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            if(i==j)con[i][j]=0;
            else if(id[i][j]==-1)dis[i][j]=inf;
            else dis[i][j]=con[i][j];
        }
    }
    for(int k=0;k<n;k++){
        for(int i=0;i<n;i++){
            for(int j=0;j<n;j++){
                dis[i][j]=max(dis[i][j],dis[i][k]+dis[k][j]);
            }
        }
    }
}
double f(int s,int t,int aim1,int aim2,double res,int flag){//p[s]到p[t]的概率,当前在s点余剩res路程可走
    if(flag==1){//未到医院
        if(con[s][t]+dis[t][aim1]+dis[aim1][aim2] > res-eps)return 0;
        else return u1*pow(eg,dis[t][aim1])+u2*pow(eg,num[s][t]);
    }
}
point choose(plane x){//给x选下一个点
    if((rand()%100)/100.0<P_jump){//突变跳
        point ret;
        ret.id=-1;
        for(int i=0;i<nump;i++){
            if(id[i][x.now.id]==-1 || i==x.now.id )continue;
            if(ret.id==-1 || con[ret.id][x.now.id]>con[i][x.now.id]){
                ret=p[i];
            }
        }
        return ret;

    }
    //---------------------------正常轮盘赌生成-------------------------
    double pob[maxn];
    pob[0]=0;
    int len=u[x.now.id].size();
    for(int i=0;i<len;i++){//生成轮盘
        point temp=p[u[x.now.id][i]];
        pob[i+1]=pob[i]+ f( x.now.id, temp.id, x.s, x.t, x.res, x.flag);
    }
    double pik=(rand()%(int)(100*pob[len]))/100.0;//转动
    for(int i=0;i<len;i++){//找轮盘停点
        if(pik<pob[i+1])return p[u[x.now.id][i]];
    }
    return p[u[x.now.id][len-1]];
}
void step(plane &x){//x向下走一步
    if(x.now==p[x.t]){//卸货
        x.flag=0;
        x.w=0;
    }
    point net=choose(x);
    //-----------------------x从now走向net
    x.way.push_back(edge(x.now.id, net.id, con[x.now.id][net.id] ) );
    num[x.now.id][net.id]++;
    x.res-=con[x.now.id][net.id];
    x.now=net;
}
double cla(){//计算当前的覆盖率
    double ret=0;
    for(int i=0;i<nume;i++){
        int ts=e[i].s,tt=e[i].t;
        if(num[ts][tt]>0&&id[ts][tt]!=-1)ret+=con[ts][tt];
    }
    return ret/sum_len;
}
void run(int x){//把x号飞机飞一遍
    if(!ant[x].way.empty()){//清空上一次飞行路径
        for(int i=0;i<ant[x].way.size();i++){
            edge temp=ant[x].way[i];
            num[temp.s][temp.t]--;
            num[temp.t][temp.s]--;
        }
        ant[x].way.clear();
    }
    ant[x].flag=1;
    do{//当x不在起点或者花费时间不足最大时间的一半时继续飞
        choose(ant[x]);
    }while(ant[x].now!=p[ant[x].s]||ant[x].cost()*2<max_cost);
}
int main(){
    init();
    input();//输入数据
    floyd();//预处理最短路
    for(int i=0;i<numa;i++){//所有飞机起飞
        run(i);
        max_cost=max(max_cost,ant[i].cost());
    }
    int T=100;
    double best=0;
    while(T--){//多次迭代维护最大覆盖率
        best=max(best,cla());
        int x=rand()%numa;
        run(x);
        printf("%lf %lf\n",cla(),best);
    }
    return 0;
}
/* 输入个数规定： n m//点数 边数 x1 y1 x2 y2 ... xn yn */