数学建模模拟退火

发布时间 2023-08-28 01:45:05作者: DLSQS
#include<bits/stdc++.h>  
using namespace std;  
#define int long long  
#define endl '\n'  
double inf = 1e8;  
double R = 6371.393;//地球半径  
double pai = 3.1415926535;  
double dd[31][31];  
int cnt = 0;  
double sum;  
struct point {  
    double lo, la;//经纬度  
    int i;  
    int kk;  
}p[31];//共有30个点  
vector<point> pp;  
double hav(double x) {//半正矢公式  
    return (double)(pow(sin(x / 2.0), 2.0));  
}  
double rad(double x) {  
    return (double)(x * pai * 1.0 / 180.0);  
}  
void dis(int x, int y) {  
    double temp = hav(p[y].la - p[x].la) + cos(p[x].la) * cos(p[y].la) * 1.0 * hav(p[y].lo - p[x].lo);  
    dd[x][y] = 2.0 * R * asin(sqrt(temp)) * 1000.0;  
    dd[y][x] = dd[x][y];  
}  
double get_sum(vector<point> p) {//更新和  
    double sum = 0;  
    for (int i = 2; i < p.size(); i++)sum += dd[p[i].i][p[i - 1].i];  
    sum += dd[p[0].i][p[p.size()-1].i];  
    //cout << sum << endl;  
    return sum;  
}  
void monte() {  
    vector<point> cur = pp;  
    for (int i = 0; i < 8e3; i++) {  
        for (int j = 0; j < 30; j++) {  
            int k = rand() % 30;  
            swap(cur[j], cur[k]);  
        }  
        double new_sum = get_sum(cur);  
        if (new_sum < sum) {  
            sum = new_sum;  
            pp = cur;  
        }  
    }  
}  
double e = 1e-16, at = 0.9999999999999, T = 1.0;  
int L = 1e5;  
void SA() {//模拟退火算法  
    while (L--) {  
        vector<point> cur = pp;  
        int c1 = rand() % 30, c2 = rand() % 30;  
        if (c1 == c2) {  
            L++; continue;  
        }  
        swap(cur[c1], cur[c2]);  
        double df = get_sum(cur) - sum;  
        double sj = rand() % 10000;  
        sj /= 10000;  
        if (df < 0) {  
            pp = cur;  
            sum += df;  
        }  
        else if (exp(-df / T) > sj) {  
            pp = cur;  
            sum += df;  
        }  
        T *= at;  
        if (T < e)break;  
    }  
}  
void show() {//将答案输出  
    sum += dd[pp[0].kk][pp[29].kk];//首位相连  
    printf("%.6f\n", sum);  
    for (int i = 0; i < 30; i++) {  
        if (i)cout << "->";  
        cout << pp[i].kk;  
    }  
    cout << "->" << pp[0].kk;  
}  
void solve() {  
    srand((unsigned)time(NULL));  
    for (int i = 0; i < 30; i++) {  
        double x, y;  
        cin >> x >> y;  
        p[i].lo = rad(x);  
        p[i].la = rad(y);  
        point temp;  
        temp.i = i;  
        temp.kk = i;  
        //cout << i << endl;  
        pp.push_back(temp);  
    }  
    for (int i = 0; i < 30; i++) {  
        for (int j = i + 1; j < 30; j++) {  
            dis(i, j);//计算第i个点和第j个点的距离  
            //cout << dd[i][j] << " ";  
        }  
        //cout << endl;  
    }  
    sum = get_sum(pp);  
    monte();  
    SA();  
    //for (int i = 1; i <= 30; i++)cout << xx[i] << endl;  
    show();  
}  
signed main() {  
    ios::sync_with_stdio(false);  
    cin.tie(0); cout.tie(0);  
    int t = 1;  
    //cin >> t;  
    while (t--)solve();  
    return 0;  
}