import matplotlib.pyplot as plt
import numpy as np
np.random.seed(2025)

Let's start by defining a planning area with a given height and width.

# Define the area dimensions
area_width = 1200
area_height = 1200

In this area, we are asked to install some cellular base stations (BS) to give coverage to the users present in the area.

We are allowed to install BS in certain fixed points, which we call candidate sites (CS). Generally speaking, the location of these sites come from agreements between the operator and the owner of the locations. CSs can be found on top of high buildings, for instance.

In our simulations, we assume that CSs are uniformely distribuited at random in the planning area. This is done for the sake of simplicity.

num_cs = 20
cs_x_coords = np.random.uniform(0, area_width, num_cs)
cs_y_coords = np.random.uniform(0, area_height, num_cs)

Installing a BS in any of our CSs has a cost. We call this installation or activation cost. As you can imagine, this cost varies a lot, according to several factors. In our simulation, we assign a random cost to every CS.

installation_costs = np.random.uniform(1, 5, num_cs)

Let's plot our planning area with the cost of every CS

plt.figure()
plt.scatter(cs_x_coords, cs_y_coords, marker='o', color='blue', label='Candidate Sites')

for i in range(num_cs):
    plt.text(cs_x_coords[i] + 20, cs_y_coords[i] + 20, f'CS {i}, c={installation_costs[i]:.2f}', color='black')

plt.xlim(0, area_width)
plt.ylim(0, area_height)
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.title('Candidate sites poisition in the planning area')
plt.legend()
plt.grid(True)

plt.show()

Now let's talk about the users in the area. We call them test points (TP). Once again, generally speaking, the position of the TPs can be any. Usually it is influenced by user position concentration in long term, or particularly important positions where the operator is highly interested in ensuring good coverage.

Same as before, we uniformely distribute them in the planning area.

num_tp = 100
tp_x_coords = np.random.uniform(0, area_width, num_tp)
tp_y_coords = np.random.uniform(0, area_height, num_tp)

Let's plot their position together with the CS.

plt.figure()
plt.scatter(cs_x_coords, cs_y_coords, marker='o', color='blue', label='Candidate Sites')
plt.scatter(tp_x_coords, tp_y_coords, marker='x', color='red', label='Test Points')

# Set plot limits and labels
plt.xlim(0, area_width)
plt.ylim(0, area_height)
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.title('Candidate sites and test point poisition in the planning area')
plt.legend()
plt.grid(True)

# Display the plot
plt.show()

In this planning instance, our planning goal is to cover as many TPs as we can by activating BSs in CSs. But our budget is limited, so we cannot just activate all of them.

To keep our budget meaningful and feasible, we set it to the average cost. This is done only for the sake of simplicity.

budget = np.mean(installation_costs)*3

But how do we define coverage?

In this planning instance, we consider a TP as covered when its received power $P_T$ is above a certain threshold $\gamma$. Note that in this planning instance we are interested only in the downlink direction.

We want to evaluate what is the received power $P_R$ at each TP according to the following formula:

$P_R = P_T+g_T+g_R-L_P,$

As we know:

• $P_T$ is the transmit power of the base stations
• $g_T$ is the transmit gain of the base stations
• $g_R$ is the receive gain of the user's terminals
• $L_P$ is the pathloss

For the sake of simplicity, we assume both receive and transmit antennas to be isotropic, with $g_T=g_R=1$ in all directions.

We will use the Okomura-Hata model to compute the pathloss in dB:

$L_P = 69.55 + 26.26\log_{10}f - 13.82\log_{10}h_B - C_H + [44.9 - 6.55\log_{10}h_B]\log_{10} d,$

where:

• $h_B [m]$ is the height of the base station antenna
• $h_M [m]$ is the hight of the user's antenna
• $f [MHz]$ is the carrier frequency
• $d [km]$ is the distance between the base station and the user
• $C_h$ is a correction factor defined as follows, using the previous parameters

$C_h = 0.8 + (1.1\log_{10}f - 0.7)h_M - 1.56\log_{10}f$

Let's define some parameters and the pathloss function...

p_t = 35 # dBm
g_t = 1
g_r = 1
gamma = -70 # dBm

f = 700 # MHz

h_b = 30
h_m = 1.3

def pathloss_oh(h_b, h_m, f, d):
  c_h = 0.8 + (1.1*np.log10(f) - 0.7)*h_m - 1.56 * np.log10(f)
  l_p = 69.55 + 26.26 * np.log10(f) - 13.82*np.log10(h_b) - c_h + (44.9 - 6.55*np.log10(h_b))*np.log10(d)
  return l_p

Let's take this chance to see how the Okomura-Hata model looks like when it's plotted

distances = np.arange(1, 11)
pathloss_values = [pathloss_oh(h_b, h_m, f, d) for d in distances]

plt.figure()
plt.plot(distances, pathloss_values)
plt.xlabel('Distance (km)')
plt.ylabel('Pathloss (dB)')
plt.title('Pathloss vs. Distance (Okumura-Hata Model)')
plt.grid(True)
plt.show()

We can finally compute $P_T$ for all test points. We start from the CS-TP distance (converted to km), and then apply the channel model formula

p_r = np.zeros((num_cs, num_tp))
for i in range(num_cs):
    for j in range(num_tp):
        d = np.sqrt((cs_x_coords[i] - tp_x_coords[j])**2 + (cs_y_coords[i] - tp_y_coords[j])**2) / 1000
        l_p = pathloss_oh(h_b, h_m, f, d)
        p_r[i, j] = p_t + g_t + g_r - l_p

Let's plot some statistics...

avg_p_r = np.mean(p_r, axis=0)

plt.figure()
plt.hist(p_r.flatten(), bins=50, density=True, alpha=0.6, color='b', label='PDF')
plt.axvline(np.mean(p_r), color='r', linestyle='dashed', linewidth=1, label=f'Average p_r: {np.mean(p_r):.2f}')
plt.axvline(gamma, color='g', linestyle='dashed', linewidth=1, label=f'p_r = gamma: {gamma}')

plt.xlabel('Received Power (dBm)')
plt.ylabel('Probability Density')
plt.title('PDF of Received Power')
plt.legend()
plt.grid(True)
plt.show()

count = 0
for j in range(num_tp):
    for i in range(num_cs):
        if p_r[i, j] > gamma:
            count += 1
            break # Move to the next test point once a suitable cs is found

print(f"Number of TPs that can be covedered by at least 1 CS: {count} out of {num_tp}")

Number of TPs that can be covedered by at least 1 CS: 100 out of 100

Let's see a subset of possible connections...

tp_indexes = [99]

plt.figure()

plt.scatter(cs_x_coords, cs_y_coords, marker='o', color='blue', label='Candidate Sites')


for i in range(num_cs):
    plt.text(cs_x_coords[i] + 20, cs_y_coords[i] + 20, f'CS {i}, c={installation_costs[i]:.2f}', color='black')

for tp_index in tp_indexes:
    plt.scatter(tp_x_coords[tp_index], tp_y_coords[tp_index], marker='x', color='red', label=f'Test Point {tp_index}')

    for i in range(num_cs):
        plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5)
        xx = (tp_x_coords[tp_index] + cs_x_coords[i]) / 2
        yy = (tp_y_coords[tp_index] + cs_y_coords[i]) / 2
        if p_r[i, tp_index] > gamma:
            clr = "green"
        else:
            clr = "red"
        plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color=clr)
        plt.text(xx, yy, f'p_r = {p_r[i, tp_index]:.2f}', fontsize=8, ha='center', va='center')

plt.xlim(0, area_width)
plt.ylim(0, area_height)
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.title('Connections between Test Points and Candidate Sites')
plt.legend()
plt.grid(True)
plt.show()

<ipython-input-33-4ee50e091759>:22: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "k-" (-> color='k'). The keyword argument will take precedence.
  plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color=clr)

We can see how, for this TP, multiple CSs can provide coverage to this TP.

The plot it's a bit messy. Let's clear it a bit.

tp_indexes = [99]

plt.figure()

plt.scatter(cs_x_coords, cs_y_coords, marker='o', color='blue', label='Candidate Sites')


for i in range(num_cs):
    plt.text(cs_x_coords[i] + 20, cs_y_coords[i] + 20, f'CS {i}, c={installation_costs[i]:.2f}', color='black')

for i in tp_indexes:
    plt.text(tp_x_coords[i] + 20, tp_y_coords[i] + 20, f'TP {i}', color='black')

for tp_index in tp_indexes:
    plt.scatter(tp_x_coords[tp_index], tp_y_coords[tp_index], marker='x', color='red', label=f'Test Point {tp_index}')

    for i in range(num_cs):
        if p_r[i, tp_index] > gamma:
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5)
          xx = (tp_x_coords[tp_index] + cs_x_coords[i]) / 2
          yy = (tp_y_coords[tp_index] + cs_y_coords[i]) / 2
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')
          # plt.text(xx, yy, f'p_r = {p_r[i, tp_index]:.2f}', fontsize=8, ha='center', va='center')

plt.xlim(0, area_width)
plt.ylim(0, area_height)
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.title('Connections between Test Points and Candidate Sites')
plt.legend()
plt.grid(True)
plt.show()

<ipython-input-34-0767d7423a74>:22: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "k-" (-> color='k'). The keyword argument will take precedence.
  plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')

To remain in budget, we would be better off selecting CS 7.

But if we add another test point that choice might not be optimal anymore...

tp_indexes = [99,90]

plt.figure()

plt.scatter(cs_x_coords, cs_y_coords, marker='o', color='blue', label='Candidate Sites')


for i in range(num_cs):
    plt.text(cs_x_coords[i] + 20, cs_y_coords[i] + 20, f'CS {i}, c={installation_costs[i]:.2f}', color='black')

for i in tp_indexes:
    plt.text(tp_x_coords[i] + 20, tp_y_coords[i] + 20, f'TP {i}', color='black')

for tp_index in tp_indexes:
    plt.scatter(tp_x_coords[tp_index], tp_y_coords[tp_index], marker='x', color='red', label=f'Test Point {tp_index}')

    for i in range(num_cs):
        if p_r[i, tp_index] > gamma:
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5)
          xx = (tp_x_coords[tp_index] + cs_x_coords[i]) / 2
          yy = (tp_y_coords[tp_index] + cs_y_coords[i]) / 2
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')
          # plt.text(xx, yy, f'p_r = {p_r[i, tp_index]:.2f}', fontsize=8, ha='center', va='center')

plt.xlim(0, area_width)
plt.ylim(0, area_height)
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.title('Connections between Test Points and Candidate Sites')
plt.legend()
plt.grid(True)
plt.show()

<ipython-input-35-be04236d9e1c>:22: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "k-" (-> color='k'). The keyword argument will take precedence.
  plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')

If we keep CS 7 active, then we need to activate another CS if we want to cover TP 90 (CS 6 is the least expensive). The toal cost would be 4.47.

There's an alternative choice: activating a subotimal CS when both TPs are taken separately, but an optimal one if they are considered together: activating only CS 10. This CS can cover both TPs with a cost of 3.59.

But should we do this? If we can't afford 4.47 (that is, CS 7 and 6), then paying 3.59 (only CS 10) is a forced choice. But if we have enough budget, keeping 2 CSs active might increase the overall coverage considering other test points...

tp_indexes = [99,90]
cs_indexes = [10]

plt.figure()

plt.scatter(cs_x_coords[cs_indexes], cs_y_coords[cs_indexes], marker='o', color='blue', label='Candidate Sites')

for i in range(num_cs):
    if i in cs_indexes:
      plt.text(cs_x_coords[i] + 20, cs_y_coords[i] + 20, f'CS {i}, c={installation_costs[i]:.2f}', color='black')

for i in tp_indexes:
    plt.text(tp_x_coords[i] + 20, tp_y_coords[i] + 20, f'TP {i}', color='black')

covered_tps = 0
for tp_index in range(num_tp):
    plt.scatter(tp_x_coords[tp_index], tp_y_coords[tp_index], marker='x', color='red', label=f'Test Point {tp_index}')

    for i in range(num_cs):
        if i not in cs_indexes:
          continue
        if p_r[i, tp_index] > gamma:
          covered_tps += 1
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5)
          xx = (tp_x_coords[tp_index] + cs_x_coords[i]) / 2
          yy = (tp_y_coords[tp_index] + cs_y_coords[i]) / 2
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')
          # plt.text(xx, yy, f'p_r = {p_r[i, tp_index]:.2f}', fontsize=8, ha='center', va='center')



plt.xlim(0, area_width)
plt.ylim(0, area_height)
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.title('Connections between Test Points and Candidate Sites')
#plt.legend()
plt.grid(True)
plt.show()
print(f'Covered TPs: {covered_tps}')

<ipython-input-36-9408886e9495>:27: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "k-" (-> color='k'). The keyword argument will take precedence.
  plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')

Covered TPs: 23

tp_indexes = [99,90]
cs_indexes = [7,6]

plt.figure()

plt.scatter(cs_x_coords[cs_indexes], cs_y_coords[cs_indexes], marker='o', color='blue', label='Candidate Sites')

for i in range(num_cs):
    if i in cs_indexes:
      plt.text(cs_x_coords[i] + 20, cs_y_coords[i] + 20, f'CS {i}, c={installation_costs[i]:.2f}', color='black')

for i in tp_indexes:
    plt.text(tp_x_coords[i] + 20, tp_y_coords[i] + 20, f'TP {i}', color='black')

covered_tps = 0
for tp_index in range(num_tp):
    plt.scatter(tp_x_coords[tp_index], tp_y_coords[tp_index], marker='x', color='red', label=f'Test Point {tp_index}')

    for i in range(num_cs):
        if i not in cs_indexes:
          continue
        if p_r[i, tp_index] > gamma:
          covered_tps += 1
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5)
          xx = (tp_x_coords[tp_index] + cs_x_coords[i]) / 2
          yy = (tp_y_coords[tp_index] + cs_y_coords[i]) / 2
          plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')
          # plt.text(xx, yy, f'p_r = {p_r[i, tp_index]:.2f}', fontsize=8, ha='center', va='center')



plt.xlim(0, area_width)
plt.ylim(0, area_height)
plt.xlabel('X-coordinate (m)')
plt.ylabel('Y-coordinate (m)')
plt.title('Connections between Test Points and Candidate Sites')
#plt.legend()
plt.grid(True)
plt.show()
print(f'Covered TPs: {covered_tps}')

<ipython-input-37-6daf040145f6>:27: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string "k-" (-> color='k'). The keyword argument will take precedence.
  plt.plot([tp_x_coords[tp_index], cs_x_coords[i]], [tp_y_coords[tp_index], cs_y_coords[i]], 'k-', alpha=0.5, color='green')

Covered TPs: 45

What is the best subset of all CSs? In other words, what is the subset of CSs that covers the most TPs with total cost below budget?

There is an optimal solution for this, but it is hard to find it.

Nontheless, the resulting network might perform pretty bad. This is because the model is too simple. It does not take into account those factors which make a good-performing mobile radio network. For instance, too many connected devices mean that the final throughput is low.

In this couse you learn how a MRN actually work, such that you can express the details of the technology into an accurate model.

And then, if you think you like working on this stuff, you can do the planning project and actually optimize the network layout, using MILP.

But we can't leave without having solved this planning instance. While getting the optimum is too much for now, let's use a simple algorithm (an heuristic) to at least get a feasible solution.

# sort cs by number of covered tps
cs_covered_tps = np.zeros(num_cs)
for i in range(num_cs):
  for tp_index in range(num_tp):
    if p_r[i, tp_index] > gamma:
      cs_covered_tps[i] += 1
cs_indices_sorted_by_covered_tps = np.flip(np.argsort(cs_covered_tps))

active_cs = []
network_cost = 0
# add the next highest-covering CS, until we run out of budget
for i in cs_indices_sorted_by_covered_tps:
  print(f'CS {i} can cover {cs_covered_tps[i]} test points')
  if network_cost + installation_costs[i] <= budget:
    active_cs.append(i)
    network_cost += installation_costs[i]
    print(f'We can add it, network cost so far: {network_cost}')
  else:
    print('We cannot add it, insufficient budget. Stopping')
    print("")
    break
  print("")

tp_is_covered = np.zeros(num_tp)
for tp_index in range(num_tp):
  for i in active_cs:
    if p_r[i, tp_index] > gamma:
      tp_is_covered[tp_index] = 1
      break
n_covered_tps = np.sum(tp_is_covered)
print(f'Covered TPs: {n_covered_tps}')

print(f'Active CS set: {active_cs}')
print(f'Final network cost: {network_cost}')
print(f'Using {(100*network_cost/budget):.2f}% of the budget')

CS 0 can cover 27.0 test points
We can add it, network cost so far: 4.7163507087626115

CS 19 can cover 26.0 test points
We can add it, network cost so far: 7.49254380468413

CS 16 can cover 25.0 test points
We cannot add it, insufficient budget. Stopping

Covered TPs: 45.0
Active CS set: [0, 19]
Final network cost: 7.49254380468413
Using 76.94% of the budget