Example 4: Network optimization

This example examines a network of multiserver queues, as shown in the figure below. Entities first arrive at node 0. After completing service at node 0, they proceed to node 1. After completing service at node 1, they proceed to node 2. After completing service at node 2, they are routed back to node 0 with probability 0.1 or exit the system with probability 0.9.

../_images/example-4.svg

The other model parameters are:

  • Five servers at each node

  • External arrivals pmf to node 0: {2: 0.3, 3: 0.7}

  • Service duration pmf at node 0: {1: 0.6, 2: 0.3, 3: 0.1}

  • Service duration pmf at node 1: {2: 0.2, 3: 0.8}

  • Service duration pmf at node 2: {1: 0.4, 2: 0.1, 3: 0.5}

  • 25 time periods

The current allocation of servers will not prevent entities from building up over time. Here's why. For steady-state analysis to be feasible, the effective offered load at each node must be less than the number of servers. The effective offered loads at nodes 1 and 2 are 8.4 and 6.3, respectively, and exceed the number of servers of five. Having a fixed time horizon is quite different from considering a steady-state scenario. For QPLEX, the effective offered load is not relevant as it can handle transient overload with ease.

To reduce the inevitable build-up of entities, suppose it is possible to add up to five servers to the network. Each server costs 1 unit, and the holding cost for each entity waiting in the buffer is 0.1 unit per period. The natural question here is:

Where to place the additional servers?

The first step is to import the correct engine and to call its constructor. The constructor is required and its argument is the number of nodes in the network. Here is the code:

from qplex import StandardNetworkMultiserver
engine = StandardNetworkMultiserver(3)

The StandardNetworkMultiserver attributes and functions are different from those of the StandardMultiserver engine. Here is the code to set its initial attributes for this example:

engine.set_number_of_servers_at_node(0, 5)
engine.set_number_of_servers_at_node(1, 5)
engine.set_number_of_servers_at_node(2, 5)

engine.set_routing_probability(0, 1, 1.0)
engine.set_routing_probability(1, 2, 1.0)
engine.set_routing_probability(2, 0, 0.1)

engine.set_number_of_external_arrivals_at_node_pmf(0, {2: 0.3, 3: 0.7})

engine.set_service_duration_at_node_pmf(0, {1: 0.6, 2: 0.3, 3: 0.1})
engine.set_service_duration_at_node_pmf(1, {2: 0.2, 3: 0.8})
engine.set_service_duration_at_node_pmf(2, {1: 0.4, 2: 0.1, 3: 0.5})

There are 56 possible incremental allocations of up to five servers to the three nodes. The maximum number of servers in the system is 20. Here is the code to execute a brute-force enumeration:

engine.mark()

server_cost = 1
waiting_cost = 0.1
initial_numbers_of_servers = [5, 5, 5]
additional_number_of_servers = 5

min_objective = float('inf')

for n0 in range(additional_number_of_servers + 1):
    for n1 in range(additional_number_of_servers - n0 + 1):
        for n2 in range(additional_number_of_servers - n0 - n1 + 1):
            incremental_number_of_servers = [n0, n1, n2]
            numbers_of_servers = \
                [initial_numbers_of_servers[i] + incremental_number_of_servers[i]
                 for i in range(engine.get_number_of_nodes())]
            objective_value = objective_function(numbers_of_servers)
            print(numbers_of_servers, objective_value)
            if (objective_value < min_objective):
                min_objective = objective_value
                minimizer = numbers_of_servers

for node in range(engine.get_number_of_nodes()):
    print('optimal number of servers at node '+str(node)+':', minimizer[node])

This code uses objective_function to calculate the cost of adding servers and the waiting costs over the 25 periods. It begins by calculating the cost of adding servers for the given allocation. It then uses the engine to obtain the pmfs of the number of entities at each node at each time epoch, from which it calculates the waiting costs over the 25 time periods. Here is its code:

def objective_function(numbers_of_servers):
    engine.restore()
    objective = 0

    for node in range(engine.get_number_of_nodes()):
        engine.set_number_of_servers_at_node(node, numbers_of_servers[node])
        objective += server_cost * numbers_of_servers[node]

    for t in range(25):
        engine.step()
        for node in range(engine.get_number_of_nodes()):
            pmf = engine.get_number_of_entities_at_node_pmf(node)
            objective += waiting_costs_single_node(pmf, numbers_of_servers[node])

    return objective

This function uses waiting_costs_single_node to calculate the waiting costs. Here is its code:

def waiting_costs_single_node(pmf, number_of_servers):
    expected_number_waiting = sum([(z - number_of_servers) * pmf[z]
                                   for z in pmf if z > number_of_servers])
    return waiting_cost * expected_number_waiting

Here is the complete code:

from qplex import StandardNetworkMultiserver

engine = StandardNetworkMultiserver(3)

engine.set_number_of_servers_at_node(0, 5)
engine.set_number_of_servers_at_node(1, 5)
engine.set_number_of_servers_at_node(2, 5)

engine.set_routing_probability(0, 1, 1.0)
engine.set_routing_probability(1, 2, 1.0)
engine.set_routing_probability(2, 0, 0.1)

engine.set_number_of_external_arrivals_at_node_pmf(0, {2: 0.3, 3: 0.7})

engine.set_service_duration_at_node_pmf(0, {1: 0.6, 2: 0.3, 3: 0.1})
engine.set_service_duration_at_node_pmf(1, {2: 0.2, 3: 0.8})
engine.set_service_duration_at_node_pmf(2, {1: 0.4, 2: 0.1, 3: 0.5})

def waiting_costs_single_node(pmf, number_of_servers):
    expected_number_waiting = sum([(z - number_of_servers) * pmf[z]
                                   for z in pmf if z > number_of_servers])
    return waiting_cost * expected_number_waiting

def objective_function(numbers_of_servers):
    engine.restore()
    objective = 0

    for node in range(engine.get_number_of_nodes()):
        engine.set_number_of_servers_at_node(node, numbers_of_servers[node])
        objective += server_cost * numbers_of_servers[node]

    for t in range(25):
        engine.step()

        for node in range(engine.get_number_of_nodes()):
            pmf = engine.get_number_of_entities_at_node_pmf(node)
            objective += waiting_costs_single_node(pmf, numbers_of_servers[node])
    return objective

engine.mark()

server_cost = 1
waiting_cost = 0.1
initial_numbers_of_servers = [5, 5, 5]
additional_number_of_servers = 5

min_objective = float('inf')

for n0 in range(additional_number_of_servers + 1):
    for n1 in range(additional_number_of_servers - n0 + 1):
        for n2 in range(additional_number_of_servers - n0 - n1 + 1):
            incremental_number_of_servers = [n0, n1, n2]
            numbers_of_servers = \
                [initial_numbers_of_servers[i] + incremental_number_of_servers[i]
                 for i in range(engine.get_number_of_nodes())]
            objective_value = objective_function(numbers_of_servers)
            if (objective_value < min_objective):
                min_objective = objective_value
                minimizer = numbers_of_servers

for node in range(engine.get_number_of_nodes()):
    print('optimal number of servers at node '+str(node)+':', minimizer[node])

And here is the output:

optimal number of servers at node 0: 5
optimal number of servers at node 1: 8
optimal number of servers at node 2: 7

This script runs in 0.732 seconds on an Apple M1 chip.