2011 MCM-B 中继站分布问题 粒子群算法实现

这次是2022美赛前的一次模拟赛，指导老师给定了2011MCM-B这道题，此题是如何分布中继站最优问题，题目先是给出了一种中继站分配带宽的系统（CTCSS），可以将带宽分为很多个 private line(PL)，通过这个可以估计出在平均人口密度下，每个中继站的覆盖半径，然后求在 $40$ 英里范围中，有 $1000$ 和 $10000$ 个用户时，如何设置中继站的位置，使得用最少的中继站将所有用户覆盖。最后讨论下在山区环境下中继站分布要注意的因素。

之前学习的粒子群算法一般用于寻找高维空间下的最值问题，每个粒子运动速度每次都带有随机性，但我们可以利用它的思路，将每一个粒子视为一个中继站，然后先均匀随机生成中继站的位置，每个粒子都向还没有被信号覆盖的用户的方向前进，距离近的移动幅度较大，距离远的移动幅度较小，因此可以定义出每个粒子的速度，再经过迭代就可以获得非常好的解。

题目

PROBLEM B: Repeater Coordination

• The VHF radio spectrum involves line-of-sight transmission and reception. This limitation can be overcome by “repeaters,” which pick up weak signals, amplify them, and retransmit them on a different frequency. Thus, using a repeater, low-power users (such as mobile stations) can communicate with one another in situations where direct user-to-user contact would not be possible. However, repeaters can interfere with one another unless they are far enough apart or transmit on sufficiently separated frequencies.
• In addition to geographical separation, the “continuous tone-coded squelch system” (CTCSS), sometimes nicknamed “private line” (PL), technology can be used to mitigate interference problems. This system associates to each repeater a separate subaudible tone that is transmitted by all users who wish to communicate through that repeater. The repeater responds only to received signals with its specific PL tone. With this system, two nearby repeaters can share the same frequency pair (for receive and transmit); so more repeaters (and hence more users) can be accommodated in a particular area.
• For a circular flat area of radius 40 miles radius, determine the minimum number of repeaters necessary to accommodate 1,000 simultaneous users. Assume that the spectrum available is 145 to 148 MHz, the transmitter frequency in a repeater is either 600 kHz above or 600 kHz below the receiver frequency, and there are 54 different PL tones available.
• How does your solution change if there are 10,000 users?
• Discuss the case where there might be defects in line-of-sight propagation caused by mountainous areas.

以下的讨论都是没有障碍物的情况

中继站半径计算

• $R$：整个用户分布半径，此题取 $40$。
• $r$：每个中继站的覆盖半径。
• $N$：用户总量。

用户密度公式：

$$\rho = \frac{N}{\pi R^2}$$

容量公式（是从另一个论文上来的）

$$\frac{3\sqrt{3}}{2}r^2\rho\leqslant 119$$

经过计算可以得出如下关系：

人口数目	1000	2000	4000	6000	8000	10000	15000
中继器半径	15.1734	10.7292	7.5867	6.1945	5.3646	4.7982	3.9178

有了半径就可以考虑覆盖算法了。

模型

蜂巢模型

这是一种非常简单的模型，可以证明用很多个相同半径的小圆覆盖大圆，完全覆盖的最优方法就是将每个小圆的圆心视为六边形的中心，六边形边长和小圆半径相同，然后将大圆放置到六边形构成的图形中，最后将六边形再换成小圆就行了。

下面这个MATLAB代码可以打出 1~13 个六边形能覆盖的最大圆。

点击显/隐代码

function main
    global R;
    global pos;
    pos = [];
    R = 40/sqrt(7);
    print(12);
    axis equal;
end

function print(n)
    if n >= 1
        drawhex(0, 0);
    end
    if n >= 2
        drawhex(sqrt(3), 0);
    end
    if n >= 3
        drawhex(sqrt(3)/2, 3/2);
    end
    if n >= 4
        drawhex(sqrt(3)/2+sqrt(3), 3/2);
    end
    if n == 5 || n == 6 || n >= 8
        drawhex(sqrt(3)+sqrt(3), 0);
    end
    if n >= 6
        drawhex(sqrt(3), 3);
    end
    if n >= 7
        drawhex(0, 3);
        drawhex(-sqrt(3)/2, 3/2);
    end
    if n >= 9
        drawhex(5*sqrt(3)/2, 3/2);
    end
    if n >= 10
        drawhex(2*sqrt(3), 3);
    end
    if n >= 11
        drawhex(3*sqrt(3)/2, 3+3/2);
    end
    if n >= 12
        drawhex(sqrt(3)/2, 3+3/2);
    end
    if n >= 13
        drawhex(5*sqrt(3)/2, 3+3/2);
    end
    if n == 1
        drawcirc(0, 0, sqrt(3) / 2);
    end
    if n == 2
        drawcirc(0, 0, sqrt(3) / 2);
        axis([-1 2.7 -1 1]);
    end
    if n == 3
        drawcirc(sqrt(3)/2, 1/2, 1);
    end
    if n == 4
        drawcirc(3*sqrt(3)/4, 3/4, sqrt(7)/2);
    end
    if n == 5
        drawcirc(sqrt(3), 3/5, 7/5);
        axis([-1 4.5 -1 1]);
    end
    if n == 6
        drawcirc(sqrt(3), 1, sqrt(3));
    end
    if n == 7 || n == 8 || n == 9
        drawcirc(sqrt(3)/2, 3/2, 2);
        if n == 9
            axis([-1.9 5.35 -1 4]);
        end
    end
    if n == 10 || n == 11
        drawcirc(sqrt(3), 3/2, sqrt(19)/2);
        if n == 10
            axis([-1.9 5.35 -1 4]);
        else
            axis([-1.9 5.35 -1 5.5]);
        end
    end
    if n == 12 || n == 13
        drawcirc(sqrt(3), 2, sqrt(7))
    end
end

function drawhex(dx, dy) 
    dx = dx * 40 / sqrt(7);
    dy = dy * 40 / sqrt(7);
    global R;
    global pos;
    pos = [pos; dx dy];
    t = 0 : pi / 100 : 2 * pi;
    x = R * sin(t);
    y = R * cos(t);
    plot(x + dx, y + dy, 'b');
    hold on
end

function drawcirc(dx, dy, r)
    dx = dx * 40 / sqrt(7);
    dy = dy * 40 / sqrt(7);
    r = r * 40 / sqrt(7);
    t = 0 : pi / 200 : 2 * pi;
    x = r * sin(t);
    y = r * cos(t);
    plot(x + dx, y + dy, 'r');
    hold on
end

下面这个代码能够打印以一个六边形为中心，画出对应的圈数下的图像

点击显/隐代码

function main
    global R;
    global pos;
    global mx;
    pos = [];
    mx = [];
    R = 5;
    dfs(6, 0, 0); // 这里就是画6圈
    t = 0 : pi / 200 : 2 * pi;
    x = 40 * sin(t);
    y = 40 * cos(t);
    plot(x, y, 'r');
    axis equal;
end

function draw(dx, dy) 
    global R;
    t = 0 : pi / 3 : 2 * pi;
    x = R * sin(t);
    y = R * cos(t);
    plot(x + dx, y + dy, 'b');
    hold on
end

function dfs(N, dx, dy)
    if N == 0
        return
    end
    global R;
    global pos;
    global mx;
    [n, ~] = size(pos);
    fg = 0;
    for i = 1 : n
        if (abs(pos(i, 1) - dx) < eps) && (abs(pos(i, 2) - dy) < eps)
            fg = 1;
            if (N <= mx(i))
                return
            else
                mx(i) = N;
            end
        end
    end
    if fg == 0
        pos = [pos; dx dy];
        mx = [mx; N];
    end
    draw(dx, dy);
    xx = R * sqrt(3) / 2;
    yy = R * 3 / 2;
    dfs(N-1, dx + xx, dy + yy);
    dfs(N-1, dx - xx, dy + yy);
    dfs(N-1, dx - xx, dy - yy);
    dfs(N-1, dx + xx, dy - yy);
    xx = sqrt(3) * R;
    dfs(N-1, dx + xx, dy);
    dfs(N-1, dx - xx, dy);
end

通过上图可以看出，要用半径为 $15.1734$ 的小圆完全覆盖 $40$ 半径的大圆，最少要用 $12$ 个。

但是现实情况下，并不一定要求全覆盖，而且因用户分布的不同，中继站的密度也会不同，所以直接用蜂巢模型并不能达到最优解。

粒子群模型

粒子群思路很简单：将每个中继站视为一个粒子，每个用户视为一个点，每次迭代时，粒子会有一个新的速度方向，每次会向这个方向移动一次。

核心：粒子只向还没有被覆盖的点方向移动，距离粒子越近的点对粒子吸引力越强，反之则越弱。（此处满足反比关系，所以公式中是利用反比例函数完成的）

• $M$：中继器的总数。
• $v_i(k)$：第 $k$ 次迭代时的第 $i$ 个中继器的运动速度。
• $x_i(k)$：第 $k$ 次迭代后的第 $i$ 个中继器的位置。
• $x_j$：第 $j$ 个发送信号的人的位置。

P.S. $x_i(k)$ 和 $x_j$ 的下标都是从 $1$ 开始的，即 $i\in[1,M]$，$j\in[1,N]$。

速度变化公式：

$$v_i(k) = w\cdot v_i(k-1)+\sum_{j\in J}a_j\frac{x_i(k-1)-x_j}{|x_i(k-1)-x_j|}$$

其中：$J$ 为当前第 $k-1$ 次迭代时，所有未被中继站覆盖到的人的位置，即

$$J=\{j\in[1,N]:|x_j-x_i(k-1)|>r, \forall i\in[1,M]\}$$

$w$ 称为惯性权重，$a_j$ 称为加速度常数，分别取值为（此处利用了反比例关系，要先减掉一个常数，使得变化较大一些）

$$w=0.6，a_j=\frac{1}{|x_i(k-1)-x_j|-r/2}$$

位置变化公式：

$$x_i(k) = x_i(k-1)+v_i(k)$$

最后为了有不同的效果还制作了两种用户生成分布的方法，一种是均匀分布，另一种是集群分布（模拟城市和乡村人口）

还加入了颜色区别，两个相交的小圆所用的颜色不能相同（避免信号冲突）

最后发现，PSO算出的结果要好很多，最下面两排为完全覆盖所需的中继站个数关系如下表

人口数目	1000	2000	4000	6000	8000	10000	15000
城市数目	4	6	7	7	7	8	8
中继站半径	15.1734	10.7292	7.5867	6.1945	5.3646	4.7982	3.9178
均匀分布	9	20	39	61	81	100	153
集群分布	8	18	35	55	74	91	140

效果图：

PSO算法部分均由C++实现：

点击显/隐代码

// calculate coverage rate
// dye each circle
// add clumped distribution
// Official version
#include <bits/stdc++.h>
#define db double
#define ll long long
#define vi vector<int>
#define pb push_back
using namespace std;
const db PI = M_PI;
const db eps = 1e-8;
const db r = 4.7982; // Radius of repeater
const db R = 40; // Radius of the whole area
const db W = 0.6; // Inertia weights
int N; // Number of repeaters
const int Nmin = 80; // change N in range [Nmin, Nmax]
const int Nmax = 100;
const int M = 10000; // Number of senders
const int T = 3000; // recursion
const int city = 8; // Number of cities
const int cR = 10; // Radius of each city
const db alpha = 0.8; // Proportion of urban population
const int cNum = M * alpha / city; // population of each city
struct Node {
	db x, y; // position
	Node(db x = 0, db y = 0):x(x), y(y) {}
	Node operator * (db c) { return Node(x * c, y * c); }
	Node operator / (db c) { return Node(x / c, y / c); }
	Node operator + (Node b) { return Node(x + b.x, y + b.y); }
	Node operator - (Node b) { return Node(x - b.x, y - b.y); }
	db abs() { return sqrt(x * x + y * y); }
} send[M]; // repeaters, senders, velocity of repeaters, best position before
vi col;
vector<Node> rep, vel;
db dis(Node a, Node b) { return sqrt((a.x-b.x) * (a.x-b.x) + (a.y-b.y)*(a.y-b.y)); }
// return random number between [0,1]
db getrand() { return 1.0 * rand() / RAND_MAX; } 
void initSend(int mode) {
	// initialize senders' position
	if (freopen("sender.txt", "r", stdin) == NULL) {
		// mode = 0 is evenly distribution
		if (!mode) {
			printf("initial evenly distribution\n");
			for (int i = 0; i < M; i++) {
				// argument (the Angle with the x-axis)
				db t = getrand() * 2 * PI;
				// magnitude (the Distance in (0, R))
				db a = sqrt(getrand()) * R;
				send[i].x = cos(t) * a;
				send[i].y = sin(t) * a;
			}
		} else {
			// mode = 1 is clumped distribution
			printf("initial clumped distribution\n");
			vector<Node> City; // save the position of cities
			for (int i = 0; i < (int)(M * alpha); i += cNum) {
				// argument (the Angle with the x-axis)
				db t = getrand() * 2 * PI;
				// magnitude (the Distance in (0, R-cR))
				db a = sqrt(getrand()) * (R - cR);
				// the center coordinate of the city
				db cx = cos(t) * a, cy = sin(t) * a;
				City.pb(Node(cx, cy));
				// Generate city citizens
				for (int j = i; j < i + cNum; j++) {
					db t = getrand() * 2 * PI;
					db a = getrand() * cR;
					send[j].x = cos(t) * a + cx;
					send[j].y = sin(t) * a + cy;
				}
			}
			for (int i = M * alpha, bool fg = 1; i < M; i++) {
				while (fg) {
					fg = 0;
					// Generate countryside citizens
					db t = getrand() * 2 * PI;
					db a = sqrt(getrand()) * R;
					send[i].x = cos(t) * a;
					send[i].y = sin(t) * a;
					// not in any city
					for (auto c : City) if (dis(c, send[i]) < cR) { fg = 1;	break; }
				}
			}
		}
	} else {
		// initialize from the file
		int m;
		scanf("%d", &m);
		assert(m == M);
		for (int i = 0; i < M; i++) scanf("%lf %lf", &send[i].x, &send[i].y);
		fclose(stdin);
	}
}
void initRep() {
	rep = vel = vector<Node>(N);
	// initialize repeaters' position evenly
	for (int i = 0; i < N; i++) {
		db t = getrand() * 2 * PI;
		db a = sqrt(getrand()) * R;
		rep[i].x = cos(t) * a;
		rep[i].y = sin(t) * a;
	}
}
void saveSender() { // save the position of senders
	freopen("sender.txt", "w", stdout);
	printf("%d\n", M);
	for (int i = 0; i < M; i++)	printf("%.3f %.3f\n", send[i].x, send[i].y);
	fclose(stdout);
}
void dye() { // dye the repeaters
	col = vi(N); // save the color id for each repeater
	for (int i = 0; i < N; i++) {
		for (int c = 1; c <= 19; c++) { // color id
			bool fg = 1;
			for (int j = 0; j < N; j++)
				if (dis(rep[i], rep[j]) < 2 * r && col[j] == c) { fg = 0; break; }
			if (fg) { col[i] = c; break; }
		}
	}
}
void saveRepeater(int num) { // save the position of repeaters
	dye();
	// Example:
	// "repeater9_1000.txt" means there are 9 repeaters and iterate 1000 times
	string s = "repeater" + to_string(N) + "_" + to_string(num) + ".txt";
	freopen(s.c_str(), "w", stdout);
	printf("%d\n", N);
	for (int i = 0; i < N; i++)
		printf("%.3f %.3f %d\n", rep[i].x, rep[i].y, col[i]);
	fclose(stdout);
}
bool cov(Node sender) { // judge whether sender has been covered
	for (int i = 0; i < N; i++)
		if (dis(sender, rep[i]) < r)
			return 1;
	return 0;
}
void PSO() {
	saveRepeater(0); // save start status
	for (int _i = 0; _i < T; _i++) {
		vector<Node> out; // senders who haven't been covered
		for (int j = 0; j < M; j++)
			if (!cov(send[j]))
				out.pb(send[j]);
		if (out.size()) {
			for (int i = 0; i < N; i++) {
				vel[i] = vel[i] * W; // multiply the weight
				Node dv;
				for (Node sender : out) {
					Node t = sender - rep[i];
					int len = t.abs();
					dv = dv + (t / (len * (len - r / 2)));
				}
				dv = dv / out.size(); // get the average
				vel[i] = vel[i] + dv; // update the velocity
				rep[i] = rep[i] + vel[i]; // move the repeater
			}
		}
	}
	saveRepeater(T); // save end status
}
db calcRate() { // calculate the coverage rate
	int tot = 0; // total covered senders
	for (int j = 0; j < M; j++)
		if (cov(send[j]))
			tot++;
	return 1.0 * tot / M;
}
signed main() {
	srand(time(NULL));
	initSend(0);
	saveSender();
	for (N = Nmin; N <= Nmax; N++) {
		initRep();
		PSO();
		int mxCol = 0;
		for (int i = 0; i < N; i++) mxCol = max(mxCol, col[i]);
		freopen("/dev/tty", "w", stdout);
		// print the number of repeaters, coverage rate and the number of colors used
		printf("%d %f %d\n", N, calcRate(), mxCol);
	}
	return 0;
}

绘图部分由MATLAB实现，只需将 sender.txt 和 repeater* 这些文件拷贝到MATLAB目录下即可，可以选择打印对应的圈的个数和迭代次数，只需选取对应的文件即可。

点击显/隐代码

% draw different colors for each circle
function main
    col = ["#F44336";"#9C27B0";"#3F51B5";"#03A9F4";"#009688";"#8BC34A";"#FFEB3B";"#FF9800";"#795548";"#607D8B";"#D81B60";"#5E35B1";"#1E88E5";"#00ACC1";"#43A047";"#C0CA33";"#FFB300";"#F4511E";"#757575"];
    R = 40;
    r = 15.17;
    %r = 4.7982;
    id = fopen('sender.txt', 'r');
    m = fscanf(id, '%d', 1);
    w = fscanf(id, '%f', [2, m]);
    x = w(1, :);
    y = w(2, :);
    plot(x, y, '.', 'color', '#546E7A');
    hold on
    t = 0 : pi / 100 : 2 * pi;
    x = cos(t);
    y = sin(t);
    plot(x * R, y * R, 'k');
    hold on
    axis equal
    for i = [3000]
		% "repeater9_1000.txt" means there are 9 repeaters and iterate 1000 times
        s = "repeater8_" + string(i) + ".txt";
        id = fopen(s, 'r');
        n = fscanf(id, '%d', 1);
        w = fscanf(id, '%f', [3, m]);
        dx = w(1, :);
        dy = w(2, :);
        c = w(3, :);
        for j = 1 : n
            plot(x * r + dx(j), y * r + dy(j), 'color', col(c(j)), 'linewidth', 1.5);
        end
    end
end

最麻烦的就是粒子群算法中速度定义，这个参数真滴难调，没调出来时候都不晓得能不能用（有点自闭了🤢），反复尝试调出来后，才发现可以的😆，而且效果不错。

最后还是要感谢队友的支持和共同努力，完成了整篇论文的写作部分😉

总结

这里总结一下这次学到的有意思的东西

两个网址：

1. AI放大图片：效果非常好，可以在放上论文前先放大一遍。

2. 好用的颜色搭配网站：在选取不同的小圈颜色时使用的（我根本不会搭配颜色，这个网站真的好用）