examples_esipp.py

# imports

# standard

import random as rand

from statistics import mean

import math

# local

import numpy as np

import networkx as nx

# imprt src.topupopt.problems.esipp as ipp

import src.topupopt.problems.esipp.utils as utils

from src.topupopt.problems.esipp.problem import InfrastructurePlanningProblem

from src.topupopt.problems.esipp.network import Arcs, Network

from src.topupopt.problems.esipp.resource import ResourcePrice

# TODO: replace this set of examples with more deterministic ones

#******************************************************************************
#******************************************************************************

def examples(solver: str,
             solver_options: dict = None,
             seed_number: int = None,
             init_aux_sets: bool = False):
    
    # test a generic mvesipp problem using the original classes
    
    # termination criteria
    
    solver_timelimit = 60
    
    solver_abs_mip_gap = 0.001
    
    solver_rel_mip_gap = 0.01

    if type(solver_options) == dict:
    
        solver_options.update({
            'time_limit':solver_timelimit,
            'relative_mip_gap':solver_rel_mip_gap,
            'absolute_mip_gap':solver_abs_mip_gap
            })
        
    else:
        
        solver_options = {
            'time_limit':solver_timelimit,
            'relative_mip_gap':solver_rel_mip_gap,
            'absolute_mip_gap':solver_abs_mip_gap
            }
    
    #**************************************************************************
    
    # no sos, regular time intervals
    
    example_generic_problem(
        solver=solver,
        solver_options=solver_options,
        use_sos_arcs=False,
        sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_NONE,
        seed_number=seed_number,
        perform_analysis=True,
        plot_results=False, # True,
        print_solver_output=False,
        irregular_time_intervals=False,
        init_aux_sets=init_aux_sets
        )
    
    # sos, cost as weight, regular time intervals
      
    example_generic_problem(
        solver=solver,
        solver_options=solver_options,
        use_sos_arcs=True,
        sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_COST,
        seed_number=seed_number,
        perform_analysis=False,
        plot_results=False,
        print_solver_output=False,
        irregular_time_intervals=False,
        init_aux_sets=init_aux_sets
        )
    
    # sos, capacity as weight, regular time intervals
                    
    example_generic_problem(
        solver=solver,
        solver_options=solver_options,
        use_sos_arcs=True,
        sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_CAP,
        seed_number=seed_number,
        perform_analysis=False,
        plot_results=False,
        print_solver_output=False,
        irregular_time_intervals=False,
        init_aux_sets=init_aux_sets
        )
    
    # sos, specific minimum cost as weight, irregular time intervals
                   
    example_generic_problem(
        solver=solver,
        solver_options=solver_options,
        use_sos_arcs=True,
        sos_weight_key=InfrastructurePlanningProblem.SOS1_ARC_WEIGHTS_SPEC_COST,
        seed_number=seed_number,
        perform_analysis=False,
        plot_results=False,
        print_solver_output=False,
        irregular_time_intervals=True,
        init_aux_sets=init_aux_sets
        )
    
    #**************************************************************************
    
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def example_generic_problem(solver: str = 'glpk',
                            solver_options: dict = None,
                            use_sos_arcs: bool = False,
                            sos_weight_key: str = 'cost',
                            seed_number: int = None,
                            perform_analysis: bool = False,
                            plot_results: bool = False,
                            print_solver_output: bool = False,
                            irregular_time_intervals: bool = False,
                            init_aux_sets: bool = False):
    
    number_periods = 2
    
    number_intraperiod_time_intervals = 5
    
    discount_rates = tuple([0.035 for p in range(number_periods)])
    
    planning_horizon = [365*24*3600 for p in range(number_periods)] # intra-period, of course
    
    if irregular_time_intervals:
        
        # TODO: adjust demand/supply levels
        
        time_step_max_relative_variation = 0.25
        
        intraperiod_time_interval_duration = [
            (planning_horizon[0]/number_intraperiod_time_intervals)*
            (1+(k/(number_intraperiod_time_intervals-1)-0.5)*
             time_step_max_relative_variation)
            for k in range(number_intraperiod_time_intervals)]
        
    else:
    
        intraperiod_time_interval_duration = [
            planning_horizon[0]/number_intraperiod_time_intervals
            for k in range(number_intraperiod_time_intervals)]
        
    # time weights
    
    # average time interval duration
    
    average_time_interval_duration = round(
        mean(
            intraperiod_time_interval_duration
            )
        )
    
    # relative weight of time period
    
    # one interval twice as long as the average is worth twice
    # one interval half as long as the average is worth half
    
    # time_weights = [
    #     [time_period_duration/average_time_interval_duration 
    #       for time_period_duration in intraperiod_time_interval_duration] 
    #     for p in range(number_periods)]
    
    time_weights = None
    
    # create problem object
    
    ipp = InfrastructurePlanningProblem(
        name='problem', 
        discount_rates={0: discount_rates}, 
        reporting_periods={0: tuple(i for i in range(number_periods))},
        time_intervals={
            0: tuple(dt for dt in intraperiod_time_interval_duration)
            },
        time_weights=time_weights
        )
    
    # add networks and systems
    
    ipp = create_generic_networks(ipp,
                                  seed_number)
    
    # set up the use of sos, if necessary
    
    if use_sos_arcs:
        
        for network_key in ipp.networks:
            
            for arc_key in ipp.networks[network_key].edges(keys=True):
                
                if ipp.networks[network_key].edges[arc_key][
                        Network.KEY_ARC_TECH].has_been_selected():
                    
                    continue
        
                ipp.use_sos1_for_arc_selection(
                    network_key, 
                    arc_key,
                    use_real_variables_if_possible=False,
                    sos1_weight_method=sos_weight_key)
    
    # instantiate
    
    ipp.instantiate(initialise_ancillary_sets=init_aux_sets)
    
    # optimise
    
    if print_solver_output:
        
        ipp.instance.pprint()
    
    out = ipp.optimise(solver_name=solver, 
                       solver_options=solver_options, 
                       output_options={},
                       print_solver_output=print_solver_output)
    
    if out:

        print('The optimisation was successful. Running post-optimisation analysis.')
        
        # run tests
        
        utils.run_mvesipp_analysis(ipp,
                                   ipp.instance,
                                   analyse_problem=perform_analysis,
                                   analyse_results=perform_analysis)
        
    else:
        
        print('The optimisation failed. Skipping results analysis.')
    
        # run tests
        
        utils.run_mvesipp_analysis(ipp,
                                   ipp.instance,
                                   analyse_problem=perform_analysis,
                                   analyse_results=False)
    
    #**************************************************************************
    #**************************************************************************
    
    # print results
    
    if plot_results:
    
        utils.plot_mves(ipp,
                        filepath='/another_folder/',
                        filename_radical='network_')
    
    #**************************************************************************
    #**************************************************************************
    
    # return something
    
    return True
                
#******************************************************************************
#******************************************************************************

def generic_problem_get_arc_techs(number_time_intervals,
                                  network_order,
                                  network_name,
                                  arc_tech_efficiencies,
                                  number_arc_technologies,
                                  peak_flow,
                                  n1,
                                  n2,
                                  distance
                                  ):
        
    min_efficiency = min(arc_tech_efficiencies.values())
    
    # note: the network order needs to be accurate
    
    capacity = [
        peak_flow*
        (1/(min_efficiency**network_order))*
        (k+1)/number_arc_technologies
        for k in range(number_arc_technologies)]
    
    min_cost = [
        (k+1)*distance*1e3*(1+rand.random())
        for k in range(number_arc_technologies)]
    
    new_arc_tech = Arcs(
        name=( 
            network_name+
            '_arc_tech_n'+str(n1)+
            '_n'+str(n2)),
        efficiency=arc_tech_efficiencies, 
        efficiency_reverse=None,
        static_loss=None,
        capacity=capacity, 
        minimum_cost=min_cost,
        specific_capacity_cost=0,
        capacity_is_instantaneous=False,
        validate=False)
    
    # return
    
    return new_arc_tech

#******************************************************************************
#******************************************************************************

def add_arc_this_way(network,
                     network_order,
                     ipp,
                     order_boost,
                     network_names,
                     g,
                     node_start,
                     node_end,
                     arc_number_key,
                     arc_tech_efficiency,
                     number_arc_technologies,
                     peak_flow,
                     distance_matrix):

    arc_tech = generic_problem_get_arc_techs(
        ipp.time_intervals[ipp.assessment_keys[0]],
        network_order[g]+order_boost,
        network_names[g],
        arc_tech_efficiency[g],
        number_arc_technologies,
        peak_flow[g],
        node_start,
        node_end,
        distance_matrix[g][node_start][node_end]
        )
                    
    # add it to the network
    
    if arc_number_key == None:
    
        network.add_directed_arc(
            node_key_a=node_start, 
            node_key_b=node_end,
            arcs=arc_tech)
        
    else:
        
        network.modify_network_arc(
            node_key_a=node_start,
            node_key_b=node_end,
            arc_key_ab=arc_number_key,
            data_dict={
                Network.KEY_ARC_TECH: arc_tech,
                Network.KEY_ARC_UND: False}
            )
    
    #**************************************************************************
    #**************************************************************************

#******************************************************************************
#******************************************************************************
        
def create_generic_networks(ipp: InfrastructurePlanningProblem,
                            seed_number: int = None):
    
    #**************************************************************************
    #**************************************************************************
    
    if seed_number == None:
        
        seed_number = rand.randint(1,int(1e5))
    
        print('Seed number: '+str(seed_number))
    
    # initialise random number generators
    
    rand.seed(a=seed_number)
    
    np.random.seed(seed=seed_number)
    
    #**************************************************************************
    #**************************************************************************
    
    # problem specification
    
    # networks
    
    min_number_networks = 2
    
    max_number_networks = 3
    
    number_networks = rand.randint(min_number_networks, max_number_networks)
    
    # network type: supply (nodes only), demand (nodes only), hybrid (both)
    
    NET_TYPE_SUPPLY = 'supply'
    NET_TYPE_DEMAND = 'demand'
    NET_TYPE_HYBRID = 'hybrid'
    
    NET_TYPES = [NET_TYPE_SUPPLY,
                 NET_TYPE_DEMAND,
                 NET_TYPE_HYBRID]

    network_type = [NET_TYPES[rand.randint(0, len(NET_TYPES)-1)]
                    for g in range(number_networks)]
    print(network_type)
    # TODO: delete print above
    # order of network
    
    min_network_order = 2 # has to be at least 2 for hybrid mode
    
    max_network_order = 4
    
    network_order = [rand.randint(min_network_order,max_network_order)
                     for g in range(number_networks)]
    
    # import and export nodes
    
    # import nodes are needed with insuf. supply
    
    min_number_import_nodes = [(1 if (network_type[g] == NET_TYPE_DEMAND or
                                      network_type[g] == NET_TYPE_HYBRID)
                                else 0)
                               for g in range(number_networks)]
        
    max_number_import_nodes = [min_number_import_nodes[g]+rand.randint(0,1)
                               for g in range(number_networks)]
    
    # export nodes are needed with insuf. demand
                               
    min_number_export_nodes = [(1 if (network_type[g] == NET_TYPE_SUPPLY or
                                      network_type[g] == NET_TYPE_HYBRID)
                                else 0)
                               for g in range(number_networks)]
    
    max_number_export_nodes = [min_number_export_nodes[g]+rand.randint(0,1)
                               for g in range(number_networks)]
    
    min_number_other_nodes = 3
    
    max_number_other_nodes = 6
    
    number_import_nodes = [rand.randint(min_number_import_nodes[g],
                                        max_number_import_nodes[g])
                           for g in range(number_networks)]
    
    number_export_nodes = [rand.randint(min_number_export_nodes[g],
                                        max_number_export_nodes[g])
                           for g in range(number_networks)]
    
    number_other_nodes = [rand.randint(min_number_other_nodes,
                                       max_number_other_nodes)
                          for g in range(number_networks)]
    
    number_nodes = [2**network_order[g]+
                    number_import_nodes[g]+ 
                    number_export_nodes[g]+
                    number_other_nodes[g]
                    for g in range(number_networks)]
    
    # arc technologies
    
    min_number_arc_technologies = [1 for g in range(number_networks)]
    
    max_number_arc_technologies = [6 for g in range(number_networks)]
    
    number_arc_technologies = [
        rand.randint(min_number_arc_technologies[g],
                     max_number_arc_technologies[g])
        for g in range(number_networks)]
    
    #**************************************************************************
    #**************************************************************************
    
    # generate data
        
    network_names = ['grid_'+str(g)
                     for g in range(number_networks)]
        
    # import prices (could be an empty dict)
    
    import_prices = {
        (g,n):[rand.random() 
               for k in ipp.time_intervals[ipp.assessment_keys[0]]]
        for g in range(number_networks)
        for n in range(number_import_nodes[g])}
    
    # export prices (lower than import ones; random if no imports prices exist)
    
    export_prices = {
        (g,n):[min(import_prices[(g,n_imp)][k]
                   for n_imp in range(number_import_nodes[g]))*
               rand.random() if number_import_nodes[g] != 0 else rand.random()
               for k in range(
                       len(ipp.time_intervals[ipp.assessment_keys[0]])
                       )]
        for g in range(number_networks)
        for n in range(number_export_nodes[g])}
    
    # static supply (negative) or demand (positive)
        
    base_flow = {
        (g,n):[rand.random() if network_type[g] == NET_TYPE_DEMAND else
               -rand.random() if network_type[g] == NET_TYPE_SUPPLY else 
               -1+2*rand.random()
               for k in ipp.time_intervals[ipp.assessment_keys[0]]]
        for g in range(number_networks)
        for n in range(number_other_nodes[g])}
    
    # positions
    
    position_nodes = [
        [(rand.random(),rand.random()) for n in range(number_nodes[g])]
        for g in range(number_networks)]
    
    # distance
    
    def distance_function(x1,x2,y1,y2):
        
        return np.sqrt((x1-x2)**2+(y1-y2)**2)
        
    distance_matrix = [
        [[distance_function(position_nodes[g][n1][0],
                            position_nodes[g][n2][0],
                            position_nodes[g][n1][1],
                            position_nodes[g][n2][1])
          for n2 in range(number_nodes[g])]
         for n1 in range(number_nodes[g])]
        for g in range(number_networks)]
    
    # determine peak demand
    
    peak_flow = [
        sum(
            max(
                [abs(base_flow[(g,n)][k]) 
                 for k in range(
                         len(ipp.time_intervals[ipp.assessment_keys[0]])
                         )
                 ]
                )    
            for n in range(number_other_nodes[g])
            )
        for g in range(number_networks)
        ]
    
    # arc tech efficiency per arc tech and grid
    
    # arc_tech_efficiency = [
    #     [1-rand.random()*rand.random()*rand.random()
    #      for k in range(ipp.number_intraperiod_time_intervals)]
    #     for g in range(number_networks)]
        
    arc_tech_efficiency = [
        {(q,k): 1-rand.random()*rand.random()*rand.random()
         for q in ipp.assessment_keys
         for k in range(ipp.number_time_intervals[q])}
        for g in range(number_networks)
        ]
                
    #**************************************************************************
    #**************************************************************************
    
    # for each network:
    # 1) create network using networkx's graph creators
    # 2) add random data to the existing nodes and edges
    # 3) add import, export and other nodes (including the relevant data)
    # 4) add arcs from these new nodes to the other nodes in the network
            
    # list of Network objects
    
    for g in range(number_networks):
        
        #**********************************************************************
        
        order_boost = (
            2 if number_import_nodes[g] and number_export_nodes[g] else 1 
            if number_import_nodes[g] or number_export_nodes[g] else 0)

        #**********************************************************************
        
        # 1) create network using networkx's graph creators
        
        if network_type[g] == NET_TYPE_DEMAND:
            
            # consumer network (positive SB_glk)
        
            new_network = Network(
                nx.binomial_tree(network_order[g],
                                 create_using=nx.MultiDiGraph))
            
        elif network_type[g] == NET_TYPE_SUPPLY:
            
            # producer network (negative SB_glk)
                
            new_network = Network(
                nx.binomial_tree(network_order[g],
                                 create_using=nx.MultiDiGraph))
            
            # reverse arc directions
            
            arc_list = []
            for arc in new_network.edges():
                arc_list.append(arc)
            for arc in arc_list:
                new_network.remove_edge(arc[0],arc[1])
                new_network.add_edge(arc[1],arc[0])
            
        else: # hybrid
        
            # join one supply grid with one demand grid
            
            G1 = nx.binomial_tree(network_order[g]-1,
                                  create_using=nx.MultiDiGraph)
                                  
            G2 = nx.binomial_tree(network_order[g]-1,
                                  create_using=nx.MultiDiGraph)
            
            nn_g1 = G1.number_of_nodes()
            arc_list = [arc for arc in G2.edges()]
            node_list = [node_key for node_key in G2.nodes()]
            for arc in arc_list:
                G2.remove_edge(arc[0],arc[1])
                G2.add_node(arc[0]+nn_g1)
                G2.add_node(arc[1]+nn_g1)
                G2.add_edge(arc[1]+nn_g1,arc[0]+nn_g1)
            for node in node_list:
                G2.remove_node(node)
                
            G = nx.union(G1,G2)                
            G.add_edge(nn_g1,0) # G2 is the supply network, G1 is the demand 1
            new_network = Network(G)
            
        
        # define the nodes as not being import, nor export nor other nodes
        
        for n in new_network.nodes:
            
            new_network.add_waypoint_node(node_key=n)
            
        # add arc data
        
        for edge in new_network.edges(keys=True):
            
            # add arc
            
            add_arc_this_way(new_network,
                             network_order,
                             ipp,
                             order_boost,
                             network_names,
                             g,
                             edge[0],
                             edge[1],
                             edge[2],
                             arc_tech_efficiency,
                             number_arc_technologies[g],
                             peak_flow,
                             distance_matrix)
            
        #**********************************************************************
        
        # compute the number of outgoing arcs per node
        
        dict_number_outgoing_arcs = {
            node: len(nx.edges(new_network,node))
            for node in new_network.nodes()}
        
        # list the nodes ordered by descending number of outgoing arcs
        
        list_nodes_descending_number_outgoing_arcs = sorted(
            dict_number_outgoing_arcs,
            key=dict_number_outgoing_arcs.get,
            reverse=True)
                        
        # list the nodes ordered by ascending number of outgoing arcs
            
        list_nodes_ascending_number_outgoing_arcs = sorted(
            dict_number_outgoing_arcs,
            key=dict_number_outgoing_arcs.get)
        
        # compute the number of incoming arcs per node
        
        dict_number_incoming_arcs = {
            node: len([node_source 
                       for node_source in new_network.predecessors(node)])
            for node in new_network.nodes()}
        
        # list of nodes ordered by descending number of incoming arcs
        
        list_nodes_descending_number_incoming_arcs = sorted(
            dict_number_incoming_arcs,
            key=dict_number_incoming_arcs.get,
            reverse=True)
                
        # list of nodes ordered by ascending number of incoming arcs
            
        list_nodes_ascending_number_incoming_arcs = sorted(
            dict_number_incoming_arcs,
            key=dict_number_incoming_arcs.get)
          
        #**********************************************************************
        
        # add import nodes
        
        for n in range(number_import_nodes[g]):
            
            # define key
            
            node_key = new_network.number_of_nodes()
            
            # res_pri = ResourcePrice(prices=import_prices[(g,n)],
            #                         volumes=None)
            
            new_network.add_import_node(
                node_key=node_key,
                prices={
                    (q,p,k): ResourcePrice(
                        prices=import_prices[(g,n)][k],
                        volumes=None
                        )
                    for q in range(ipp.number_assessments)
                    for p in range(ipp.number_reporting_periods[q])
                    for k in range(ipp.number_time_intervals[q])
                    }
                )
            
            # add arc from import node to a node with many outgoing arcs
            
            add_arc_this_way(new_network,
                             network_order,
                             ipp,
                             order_boost,
                             network_names,
                             g,
                             node_key,
                             list_nodes_descending_number_outgoing_arcs[n],
                             None,
                             arc_tech_efficiency,
                             number_arc_technologies[g],
                             peak_flow,
                             distance_matrix)
                    
        #**********************************************************************
            
        # add export nodes
        
        for n in range(number_export_nodes[g]):
            
            # define key
            
            node_key = new_network.number_of_nodes()
            
            # res_pri = ResourcePrice(prices=export_prices[(g,n)],
            #                         volumes=None)
            
            new_network.add_export_node(
                node_key=node_key,
                prices={
                    (q,p,k): ResourcePrice(
                        prices=export_prices[(g,n)][k],
                        volumes=None
                        )
                    for q in range(ipp.number_assessments)
                    for p in range(ipp.number_reporting_periods[q])
                    for k in range(ipp.number_time_intervals[q])
                    
                    }
                )
            
            # add arc from node with many incoming arcs to the export node
            
            add_arc_this_way(new_network,
                             network_order,
                             ipp,
                             order_boost,
                             network_names,
                             g,
                             list_nodes_descending_number_incoming_arcs[n],
                             node_key,
                             None,
                             arc_tech_efficiency,
                             number_arc_technologies[g],
                             peak_flow,
                             distance_matrix)
            
        #**********************************************************************
        
        # identify import and export nodes
        
        new_network.identify_node_types()
        
        #**********************************************************************
        
        # demand/supply nodes: create them and add arcs to random nodes
        
        demand_node_counter = 0
        
        supply_node_counter = 0
        
        for n in range(number_other_nodes[g]):
            
            # add demand/supply node
            
            node_key = new_network.number_of_nodes()
            
            new_network.add_source_sink_node(
                node_key=node_key,
                # base_flow=base_flow[(g,n)],
                base_flow={
                    (q, k): base_flow[(g,n)][k]
                    for q in ipp.assessment_keys
                    for k in range(
                            len(ipp.time_intervals[q])
                            )
                    }
                )
            
            # differentiate node placement based on the static flow needs
            
            # this will tend to ensure feasibility
            
            if min(base_flow[(g,n)]) >= 0:
                
                # demand node: 
                # from nodes with zero/few outgoing arcs to the demand node
                
                node_key_start = list_nodes_ascending_number_outgoing_arcs[
                    demand_node_counter]
                
                add_arc_this_way(new_network,
                                 network_order,
                                 ipp,
                                 order_boost,
                                 network_names,
                                 g,
                                 node_key_start,
                                 node_key,
                                 None,
                                 arc_tech_efficiency,
                                 number_arc_technologies[g],
                                 peak_flow,
                                 distance_matrix)
                
                # increment counter
        
                demand_node_counter = demand_node_counter + 1
                
            elif max(base_flow[(g,n)]) <= 0:
                
                # supply node: 
                # from the supply node to nodes with zero/few incoming arcs
                
                node_key_end = list_nodes_ascending_number_incoming_arcs[
                    supply_node_counter]
                
                add_arc_this_way(new_network,
                                 network_order,
                                 ipp,
                                 order_boost,
                                 network_names,
                                 g,
                                 node_key,
                                 node_key_end,
                                 None,
                                 arc_tech_efficiency,
                                 number_arc_technologies[g],
                                 peak_flow,
                                 distance_matrix)
                
                # increment counter
                
                supply_node_counter = supply_node_counter + 1
            
            else:
                
                # demand/supply node
                
                # add two arcs:
                # arc 1) from an import node or nodes directly or indirectly
                # connected to an import node (from which imports are possible)
                # arc 2) to an export node or nodes directly or indirectly co-
                # nnected to an export node (from which exports are possible)
                
                #**************************************************************
                
                # arc 1
                
                # randomly select a starting node
                
                # for each import node
                
                for import_node in new_network.import_nodes:
                    
                    # select a node with few outgoing arcs
                
                    node_key_start = list_nodes_ascending_number_outgoing_arcs[
                        supply_node_counter]
                    
                    # check if there is a path between them
                    
                    if nx.has_path(new_network, import_node, node_key_start):
                        
                        # call random for comparison purposes
                        
                        rand.randint(0,1)
                        
                        # update the counter
                        
                        supply_node_counter = supply_node_counter + 1
                        
                        # if there is, break
                        
                        break
                    
                    # if not, continue
                
                    # TODO: while loop to try more times with each import node
                
                else: 
                    
                    # randomly select an import node
                    
                    node_key_start = new_network.import_nodes[
                        rand.randint(0,len(new_network.import_nodes)-1)]
                
                # add arc
                
                add_arc_this_way(new_network,
                                 network_order,
                                 ipp,
                                 order_boost,
                                 network_names,
                                 g,
                                 node_key_start,
                                 node_key,
                                 None,
                                 arc_tech_efficiency,
                                 number_arc_technologies[g],
                                 peak_flow,
                                 distance_matrix)
                
                #**************************************************************
                
                # arc 2
                
                # randomly select an end node
                
                # for each export node
                
                for export_node in new_network.export_nodes:
                    
                    # select a node with few incoming arcs
                
                    node_key_end = list_nodes_ascending_number_incoming_arcs[
                        demand_node_counter]
                    
                    # check if there is a path between them
                    
                    if nx.has_path(new_network, node_key_end, export_node):
                        
                        # call random for comparison purposes
                        
                        rand.randint(0,1)
                        
                        # update the counter
                
                        demand_node_counter = demand_node_counter + 1
                        
                        # if there is, break
                        
                        break
                    
                    # if not, continue
                
                    # TODO: while loop to try more times with each export node
                
                else: 
                    
                    # randomly select an export node