Criticality analysis of a country’s transport network via an agent-based supply chain model

Building resilience against shocks has become a pillar of sustainability. By understanding how the different components of an economy interact in times of crisis, we can design resilience strategies that go beyond building walls or dams. We formulate an original agent-based model to explore a crucial pathway through which a disaster affects the economy: the transport–supply chain nexus. The model simulates the behaviour of firms facing transport and supply disruptions and estimates the resulting indirect losses. As an illustration, the model is used to assess the criticality of Tanzanian roads, which are vulnerable to floods. We report three main results. First, the model generates maps that identify the transport infrastructure assets that are most critical for specific supply chains: roads that are most important for food security are different from those supporting international trade, for instance. Second, economic losses from transport disruptions increase non-linearly with the duration of disruptions, highlighting the benefit from fast repairs. Third, by combining economic and transport modelling, we can consider a broader spectrum of interventions. Beyond strengthening the transport system, it is also possible to make supply chains more resilient to disruption with, for instance, sourcing decisions or inventory management. Assessing the performance of an economy in times of crisis is a difficult task. This study presents an agent-based model capturing the behaviour of firms facing transport and supply disruptions due to natural disasters, and shows possible paths to reinforce infrastructure and supply chain resilience.

B uilding resilience against disasters has become a key adaptation strategy for climate change. But our ability to adequately track the impact of disasters has remained limited. Disasters are still often estimated by summing up reconstruction costs 1,2 without considering the broader implications for the functioning of the economic system. This approach gives little insight into the real economic role of specific assets and the complex propagation of the shock through supply chains. It is likely to lead to substantial underestimates of the severity of disasters 3 .
Transport systems in particular play a key role during and after disasters because of the combination of their high vulnerability and their importance for the functioning of the rest of the economic system. Because transportation systems have become extremely efficient and reliable, economies have become increasingly dependent on them. Also, because transportation infrastructure assets are by nature localized, they are affected very heterogeneously by natural disasters and, in turn, affect the rest of the economy very heterogeneously. Two bridges may be of equal engineering quality and construction cost and may experience the same damages (defined as the repair cost), but the economic implications of closing each bridge during the reconstruction might be completely different. One bridge might be playing a vital role in the local economy (for example, the Bay Bridge in California 4 ) while the other one might be redundant, making it possible to close it without notable economic impacts.
We propose a model that improves our ability to understand and assess the economic consequences of natural disasters by analysing the transport-supply chain interactions in space and time. Our model allows us to track how the disruption of transport nodes or links perturbs the flows of goods in supply chains and how these perturbations then affect households, firms and trade. This model is a new tool for evaluating the implications of transport disruptions, identifying bottlenecks and planning for resilience. It innovates by bridging two recent developments, namely, transport network study and firm-level input-output analysis, that used to be considered in isolation.
In the last decade, increasing efforts have been made to quantify the resilience of transport systems to natural disasters and use this information to inform decision makers. However, most of the available tools to assess the impact of disruption in transport systems rest on network theory and the concept of percolation 5,6 . Based on data on vehicles, passengers and freight, impacts are measured by the increase in travel time 7 , loss of accessibility 8 or quantities of blocked flows 9,10 . Such indicators can then inform infrastructure planning 11 .
But the economic impact of disrupted transport infrastructure does not end there. Transport disruptions translate into supply chain disruptions that act as amplifiers of disaster-induced economic shocks 12 . Multiple case studies provide detailed accounts of the propagation mechanisms, including studies of the 2011 Tōhoku earthquake in Japan [13][14][15] and the 2011 floods in Thailand 16,17 . These mechanisms are often studied using input-output models [18][19][20] , computable general equilibrium models (CGEs) 21 or hybrid models 22 . Some of these models explicitly represent the transport sectors to capture the impact of transport infrastructure destruction on the rest of the economy 23,24 .
However, the level of aggregation of standard models is problematic, and the explicit representation of firm-level interactions is preferable 25 . According to case studies, the cascades of disruptions are due to specific patterns of interdependencies between firms rather than sectors 15,26,27 . Heterogeneities among firms and in production networks have the potential to amplify localized shocks [28][29][30][31][32] . Aggregated models also have limited applicability for planning purposes because they cannot capture the role of specific transport nodes or edges. This paper brings together the firm-level input-output approaches with transport network analysis. We argue that firm-to-firm models cannot satisfactorily assess the impacts of a disaster without considering how goods and services are transported. Likewise, transport system models cannot assess the criticality of a transport link without considering the type of goods and services it transports and its role in connecting suppliers and clients. Our models contribute to the literature by combining these two approaches into an integrated assessment framework.

Célian Colon 1 ✉ , Stéphane Hallegatte 2 and Julie Rozenberg 2
Building resilience against shocks has become a pillar of sustainability. By understanding how the different components of an economy interact in times of crisis, we can design resilience strategies that go beyond building walls or dams. We formulate an original agent-based model to explore a crucial pathway through which a disaster affects the economy: the transport-supply chain nexus. The model simulates the behaviour of firms facing transport and supply disruptions and estimates the resulting indirect losses. As an illustration, the model is used to assess the criticality of Tanzanian roads, which are vulnerable to floods. We report three main results. First, the model generates maps that identify the transport infrastructure assets that are most critical for specific supply chains: roads that are most important for food security are different from those supporting international trade, for instance. Second, economic losses from transport disruptions increase non-linearly with the duration of disruptions, highlighting the benefit from fast repairs. Third, by combining economic and transport modelling, we can consider a broader spectrum of interventions. Beyond strengthening the transport system, it is also possible to make supply chains more resilient to disruption with, for instance, sourcing decisions or inventory management.
Our model builds on Henriet et al. 31 . They developed a firm-to-firm model and showed that estimated losses from natural disasters were much larger when firm-level dynamics were considered. Inoue and Todo 33 extended and applied this firm-to-firm model with a comprehensive dataset of the supplier-buyer linkages of Japan. Their model better reproduces the macroeconomic trajectory that followed the 2011 Tōhoku earthquake than aggregated input-output models. Our work adds the transport network layer to this model. From a transport analysis perspective, Kim et al. 34 and Ham et al. 35 mapped flows of commodities on the highways and railways of the United States and estimated the increase in transport costs due to the disruption of highway links. Our work extends this approach by modelling the behaviour of firms and the cascades of impacts throughout supply chains. It shares similarities with the work of Tatano and Tsuchiya 36 , who analysed the impact of the disruption of the main highways and railways of Japan using a multiregional CGE model, but our model has a finer spatial and economic granularity.
The model is applied to the United Republic of Tanzania. Combining data on the transport sector from Pant et al. 9 and OpenStreetMap contributors 37 , the input-output table for the country 38 , a listing of Tanzanian firms with their basic characteristics and the results of a dedicated firm survey that collected data on inventories and vulnerability to disruptions 39 , we created a supply chain model that fully embeds the transport network 40 . A set of firmsone per sector-is modelled on each important node of the transport network, along with households, who buy firms' outputs. Each firm buys inputs from and sells outputs to other firms located in the other nodes of the transport network. Firms are also connected to international suppliers and buyers, modelled at the country's main entry and exit points. Extended Data Fig. 1 illustrates the coupling of the supply chain and transport networks.
In normal times, suppliers provide intermediate goods to other firms using the least-cost path on the road network. This journey is associated with a transport cost that is paid for by the client. When a road gets disrupted, two scenarios can unfold. If there is an alternative, more expensive path, the increase in transport cost is transferred to the client and propagated further, leading to an increase in the final price paid by consumers. If there is no other path and the client cannot receive goods from the supplier on time, then the lack of input can (1) be compensated for by other suppliers, (2) be smoothed using inventories or (3) reduce the client's production. Shortages then propagate along supply chains and may eventually lead to drops in final consumption.
The final macroeconomic cost (that is, indirect losses) is estimated as the sum of two impacts: (1) the increased cost for final consumers, which is equal to the sum of all the incremental transportation costs in the economic system, and (2) the value of 'missing goods and services' valued at their pre-shock price, which is a conservative estimate of the impact on consumers and therefore of the welfare losses. At this stage, the model does not distinguish between consumers in a given location (for example, by income class) and does not represent consumers' behavioural response when prices change or goods are not available. We use our model to assess the criticality of the road network of the United Republic of Tanzania, motivated by the vulnerability of this economy to floods caused by heavy rainfalls or tropical storms. According to the World Bank's observers, the country is subject to small, frequent floods, which often interrupt roads for a few days, but also to larger floods that severely damage the transport network, which then takes very long to get repaired. Pant et al. 9 reported 13 such sizeable transport-disrupting events between 2014 and 2016. Using a comprehensive set of scenarios, we estimate the indirect economic losses caused by road disruptions for the country's households and international buyers. We draw maps that show the criticality of each node and link in the transport network, which can help prioritize intervention to mitigate flood impacts. The main behaviour of the model is then analysed, especially the crucial role played by the duration of the transport disruptions and the differentiated vulnerability of different supply chains. The model helps assess a variety of interventions able to reduce economic losses from road disruptions and build resilience.

Results
First, we simulate a 1 week disruption of each transport infrastructure asset (that is, the nodes and edges of the transport network) one by one and show the resulting economic impacts on the transport map. When all supply chains are considered, these criticality maps show the criticality of transport infrastructure for the overall economy, including the producers and consumers located in different regions. But if we represent the losses related to selected products only, such as food, we can generate criticality maps for specific supply chains. The four maps of Fig. 1 show, for four products, the aggregate cost (for all Tanzanian consumers, irrespective of where they live) resulting from a 1 week disruption of each road and bridge. There is uncertainty in the supply chain network, which is reconstructed based on aggregate data and with some randomness in the firm-to-firm connections, so we reproduce the same exercise on ten different networks to confirm the stability of the criticality maps.
These maps highlight the first finding of this paper: the criticality of an infrastructure asset is highly dependent on the product that is considered. For example, the roads and bridges that are important for food security may not be important for services. The criticality of an asset depends on its users 3 .
For example, the criticality map of the United Republic of Tanzania for agricultural products, shown in Fig. 1a, reveals that the supply chains underpinning those products are relatively less impacted by transport disruptions. The map is different when we consider food products-here, food products consist of processed products (for example, sugar), whereas agricultural products are unprocessed (for example, sugar cane). The production and consumption of food products rely much more on the transport system. For some of the most critical roads, a 1 week disruption can induce a 1.5% to 2% decrease in national daily consumption of food products (Fig. 1b).
The northeast-southwest corridor, which connects the port of Dar es Salaam to Zambia and, further, to the Democratic Republic of Congo, is particularly critical to international trade (Fig. 1d). These roads carry large amounts of freight from and to Central African countries. Plans to redevelop (currently unreliable) rail freight between Dar es Salaam and Zambia are likely to reduce such criticality 41 .
This first result confirms the value of connecting a supply chain model to a transport model and suggests that traditional criticality analyses of the transport system, which consider the transport network but not its users, may fail to identify the right priorities for intervention.
The second general insight from these simulations is that the economic impact increases non-linearly and extremely quickly with the duration of transportation disruptions. The temporary rerouting of freight flows triggered by 1 week disruptions mostly drives indirect economic losses (Fig. 2). Thanks to inventories, there are almost no shortages or business interruptions and losses remain limited.
But if disruptions last longer, some firms start to run out of inventories and are forced to delay their deliveries, leading to cascading impacts and scarcity of final goods and services. Scarcity impacts welfare much more than an increase in prices due to rerouting. A 2 week disruption is, on average, five times costlier than a 1 week one. The longer the disruption, the further shortages propagate in the supply chains. A 4 week disruption is on average 24 times costlier than a 2 week one.
This observation leads to the third insight from this modelling exercise. Improving resilience can be done by strengthening the transport system (supply-side action), but also by making the users of this system better able to manage disruptions (demand-side action).
To illustrate these two approaches, we define three interventions for each and compare their potential for reducing the indirect losses of road disruptions (Fig. 3). We target two critical corridors identified in Fig. 1 and shown in Extended Data Fig. 2: a portion of the T1 road south of Dodoma and a portion of the T7 road on the east coast. Together, the two segments represent 181 km of trunk roads. We consider (i) strengthening them to reduce their vulnerability to floods, (ii) increasing redundancy by building new roads in the  area and (iii) improving preparedness so that road damages can be repaired more quickly, thereby reducing the duration of disruptions. Figure 3a shows that these three actions that provide similar benefits in terms of avoided final losses. Strengthening the road with wider culverts and better materials could halve the likelihood of disruption. According to Miyamoto International 42 , building a more robust road would cost about 2% to 4% of the US$1 million per km construction cost estimated for such roads 43 . Therefore, the incremental cost of the strengthening strategy would amount to US$3.6 million to US$7.2 million, provided that the strengthening is done at the construction stage (retrofitting the road would be much more expensive). Building about 200 km of new roads to offer an  alternative corridor could lead to similar benefits, at a much higher cost of about US$200 million. These two estimates are likely to be heavily dependent on the country's topography. Interestingly, the same reduction in indirect losses could be achieved without capital investment by reducing the disruption duration from 4 weeks to 3 weeks, using, for instance, better monitoring and reporting of road damages and quicker dispatch of staff and equipment. Since better monitoring and maintenance of roads is widely considered cost-effective 44 , by reducing the need for major works, this intervention may well have a negative cost even without considering resilience gains.
Another approach to reducing the cost of transport disruptions is to ensure that firms can manage them better. This resilience can be strengthened through (i) larger inventories, (ii) a larger number of suppliers per firm or (iii) reducing the distance between suppliers and buyers-that is, adopting local-sourcing strategies. Figure 3b shows the loss-mitigating potential of these measures on the two targeted corridors, while Extended Data Fig. 3 shows the impact of the same measures on the 300 most critical transport nodes. Figure 3b shows that having one extra week of inventory significantly reduces the losses triggered by the prolonged disruption of these roads, as does having two suppliers per input instead of one. Those measures are, however, most effective for long disruptions, since firms already have enough inventory to manage one-week disruptions (according to the firm survey realized for this study). It may be the case that they have such large inventories because of the regular transport disruptions in Tanzania, a response identified in other countries 45 . Extended Data Fig. 3 shows that these results are valid for all assets: holding larger inventories and having more suppliers massively reduce losses for large disasters.
In Fig. 3b, local-sourcing strategies massively reduce losses for both short and long events. When the average distance between supplier and buyer is reduced, less freight is being transported on the two considered corridors. That way, disrupting these roads generates fewer perturbations and induces fewer losses. But, as shown in Extended Data Fig. 3, this result is not valid everywhere. In general, local-sourcing strategies reduce the indirect losses of short disruptions but increase those of long ones. Relying on local suppliers and clients reduces the need for transport and exposure to distant transport-disruption disasters. But when local disaster hits all suppliers and clients simultaneously, a consequence of local sourcing, then interdependency across firms makes the recovery longer and more difficult. Had they maintained relationships with distant suppliers and clients, they would be able to source inputs and sell products as soon as transport returned to normal (something observed empirically after the 2011 earthquake in Japan 14 ).

Discussion
We combine a geospatially disaggregated firm-level input-output model with a model of the Tanzanian transport network to analyse the full economic criticality of each road segment. This coupling highlights how different supply chains have different exposures to transport disruptions and amplify those shocks differently. As a result, sector-specific criticality maps exhibit different patterns. The bottlenecks for food security will not be the same as those for trade competitiveness.
Quantifying criticality in this way can inform investment decisions. Strengthening infrastructure assets can reduce disaster losses, but with limited resources to allocate to infrastructure 46 , it is crucial to prioritize the most critical assets for the economic system's functioning. Doing so requires considering not only construction costs and transport flows but also accounting for the criticality of the assets (that is, their importance for various sectors of the economy). Overlaying hazard flood maps with our criticality maps (Fig. 1) would provide priorities for making new roads and bridges more resistant to floods.
But resilient infrastructure systems do not only rest on stronger assets. As seen in Fig. 2, each additional day of disruption gets increasingly costly. This fast and non-linear increase in indirect costs with disruption duration results from the interplay between transport continuity and inventory management. It is a general result and is not specific to the country investigated. It is consistent with Henriet et al. 31 and Inoue and Todo 33 , who show that firm-level input-output models exhibit larger cost amplification effects than sector-level ones. We obtain a similar behaviour in Extended Data Fig. 4. The fast increase also highlights how crucial inventories are in the resilience of the manufacturing sector. Building the capacity to quickly restore services on damaged roads or bridges 47 can generate the same economic benefits as building more or stronger roads-but at a lower cost.
Managing transport demand by reorganizing supply chains is also an option for coping with disasters. Such an approach is rarely considered in disaster risk management, which primarily addresses direct impacts and often focuses on the supply of infrastructure services. As expected, adding inventories and double-sourcing strategies significantly reduces indirect costs. The model can be used to identify the type of inputs for which additional safety stocks would be most beneficial for a country. Preliminary results, shown in Extended Data Fig. 5, suggest that food products are particularly critical for the United Republic of Tanzania. Such inventories are, however, costly due to storage needs and the risk of wastage, which is particularly acute for perishable goods 45 . Further analyses could help identify where such inventories are most needed and which firms are most critical in propagating impacts.
Relying on local suppliers and clients has a mixed effect. It largely reduces the need for transport and thus the exposure to disasters. But failing to keep relationships with distant partners can severely aggravate the losses. This result, which we expect not to depend on the specific country investigated, is consistent with findings related to the 2011 Tōhoku earthquake 14 and Hurricane Sandy 48 . It feeds into the globalization-versus-relocation debate, which has garnered increasing interest since the COVID-19 crisis began (see, for example, the work of Baldwin and Evenett 49 ). Our results suggest that relocating supply chains can boost resilience only if global ties are maintained to help manage local shocks.
The resilience agenda tends to focus on building more and stronger infrastructure, while many countries do not have the fiscal space to pay for all the infrastructure they need. Instead, our model highlights 'soft' interventions focused on improving the maintenance and accelerating the repairs of infrastructure assets, possibly delivering resilience at much lower costs. Furthermore, working with the users of infrastructure services, as suggested by Hallegatte et al. 3 , can also offer opportunities to build resilience at an affordable cost, even in poor environments.
Our main conclusions-the criticality maps, superlinear reaction to disruption duration and importance of comparing supply-side and demand-side interventions-are robust, but the quantified estimates from this study have large uncertainties. These uncertainties are due to the poor quality of the data, the probabilistic reconstruction of the supply chains, cut-offs in the number of firms and transport nodes and behavioural assumptions in the firm model (Supplementary Table 11). Better calibration and validation will be critical to improving the model, but this will require dedicated data collection and post-flood firm surveys to measure the firm-level impacts.
The model can also be improved in several directions. Although service sectors are included in the model, their operations are assumed to be independent of the transport network. This simplification allowed us to focus on freight but dwarfed the role played by service firms in connecting supply chains, such as third-party logistic providers. The impact of disasters on workers' ability to reach workplaces is an important source of economic losses and is not explicitly considered in the current version of the model. Finally, the model relies on simple input-output mechanisms and cost-minimizing behaviour, which do not include the (limited) possibility of substituting inputs for each other during a crisis. The importance of these mechanisms is evident when comparing different purchasing responses from firms (the 'reactivity rate' in the Supplementary Methods). And, while it was not necessary for this exploration of the current vulnerability of the United Republic of Tanzania, extending the model to transport multimodality will be crucial for its application to more complex transport systems.

Methods
We model the road network of the United Republic of Tanzania, which carries more than 99.8% of the trade flows 9 . The continental part of the country consists of 25 regions and 169 districts. We leave the island regions of Zanzibar out of the scope of the study. Spatial distributions of gross domestic product (GDP) and population are used to locate 324 major transport nodes, called origin-destination (OD) nodes 50,51 . We estimate the production per sector and per district by using a business registry built and maintained by the National Bureau of Statistics of the United Republic of Tanzania that provides data on about 257,894 companies for the year 2017, including sector, administrative district and number of employees. We retain 42 sectors, based on the Global Trade Analysis Project (GTAP) classifications 38 . Corrections are applied using country-wide data on sectoral production 38 , spatially explicit GDP data 51 and land cover data 52,53 . Final demand per district and per sector is estimated using population census and GTAP consumption data 38,54 . On each one of the 324 OD nodes, we model one representative household and, for each sector with substantial production at this node, one representative firm, leading to 1,962 representative firms. This setting allows us to represent the spatial distribution of population and activities along the national and regional roads. Because we aim at a country-level analysis, we do not model the last-mile accessibility to rural villages nor agricultural fields, which is best studied at the local level. Utilities and services are modelled differently because their outputs are not conveyed through the transport network as freight flows. We model country-wide representative firms for each one of these sectors and do not assign them to a specific OD node. Supplementary Figure 4 analyses the impact of the number of representative firms on the results.
Firms have a Leontief production function calibrated with GTAP input-output tables, as in the work of Henriet et al. 31 and Inoue and Todo 33 . As opposed to Inoue and Todo 33 , we do not have the suppliers and buyers of each firm in the data. We reconstruct the supply chain links using a modified gravity model, based on the size, sector and location. For each sector from which they need input, representative firms randomly pick s suppliers. For firm A of a relevant sector, the probability of being chosen as supplier by firm B grows with A's size and decreases with the distance between A and B. This procedure allows the construction of the firm-to-firm network (see below for a discussion of the uncertainty created by this procedure). We assume that households only buy from firms located on the same OD node. For imports and exports, we explicitly model neighbouring countries with their border crossing and groups of non-neighbouring countries. For exports, countries choose the firms from which they buy based on their size. For imports, Tanzanian firms pick countries as suppliers based on the country's volume of exports to the United Republic of Tanzania and based on their distance to the country's entry point into the United Republic of Tanzania. Last, we build the country-to-country trade links for transit flows. Each one of these supplierbuyer relationships is associated with the lowest-cost path in the road network linking the OD node of the supplier to the OD node of the buyer. Given the final demand of households and countries, we can generate the map of the traded freight flows. Supplementary Figure 5 analyses the impact of the probabilistic network reconstruction on the results and the uncertainty created by the absence of data on the real firm-to-firm network (using 60 randomly generated networks, which is more than in Fig. 1 and confirms the robustness of the results).
The model is a discrete dynamic map X t+1 = f(X t ), in which each time step represents a week. At each time step, firms plan production, order supplies, produce, adjust prices, deliver output and collect orders. These firm-level input-output dynamics build on the model of Henriet et al. 31 Firms manage their inventory to maintain a target stock, which we calibrate using data from a survey made in 2018 of 800 Tanzanian firms. In the absence of perturbations, the model is at a steady state. When a node or a link of the transport network is disrupted, we model one of two mechanisms to capture the indirect economic losses, depending on whether an alternative pathway exists.
If an alternative pathway exists, supply chain flows are rerouted through the lowest-cost journey available, leading to higher-than-usual transport costs. Suppliers transfer these higher costs onto their clients by increasing their selling price. In turn, clients transfer such an increase in input costs to their own clients. Price increases propagate down to households and international buyers, which have to spend more for the same basket of goods.
If an alternative pathway does not exist, products that were supposed to be delivered are held at the producers' premises. If prolonged beyond the flexibility offered by inventories, this situation may lead to a shortage of inputs for clients and reduced production for them. At the end of the supply chains, households and international buyers may incur consumption losses due to scarcity of some inputs (valued at pre-shock prices).
Indirect losses are the sum of two economic impacts, namely spending due to increased transport costs and consumption loss due to scarcity in some goods and services. These two impact mechanisms-more expensive deliveries accepted by clients versus no deliveries, compensated by inventories-represent the two extremes of a continuous reality, in which firms adapt their purchase volumes to prices and product availability. As in Henriet et al. 31 , we assume that supply chains are fixed in the short term, that is, that firms do not switch to new suppliers within the few weeks of the simulation; however, they can adjust the quantities ordered from existing suppliers.
Our measure of indirect losses sums both impacts into a single monetary metric. For households, these impacts do not strictly represent monetary flows, but instead represent disruption-induced inflation (and therefore reduction in purchasing power) and scarcity in some goods and services in local markets. Similarly, the indirect losses for international buyers may not represent an actual added cost if international trade is highly competitive. However, they quantify a loss of competitiveness of Tanzanian exports.
To capture uncertainties related to the reconstruction of supply chain links, we simulate disruptions on multiple supply chain networks and average the estimates. We perform criticality analyses of the transport network. They consist of disrupting each transportation infrastructure one by one and measuring the associated indirect loss. We simulate disruptions that last for 1-4 weeks. For each duration, we perform statistical analysis based on the 300 most critical infrastructures. Supplementary Figure 5 shows the dispersion of the results due to the probabilistic reconstruction of the supply chain network. It shows that the identification of the most critical transport nodes and edges is robust to variations in the reconstructed networks.
The model is fully described in the Supplementary Methods, first formulating the model and then explaining how the data were processed to calibrate it.

Data availability
Input-output tables and international trade data were extracted from the GTAP database 38 . It is a licensed database provided by GTAP with restricted shareability, so the data are not publicly available. The firm registry and population census are data provided by the authorities of the United Republic of Tanzania. They are not publicly available; requests should be addressed to the National Bureau of Statistics of the United Republic of Tanzania. Gridded population data, spatial GDP data and spatial land cover data are available online; see references for direct links in CIESIN 50 , Kummu et al. 51 and RCMRD 53 . Inventories were estimated from a World Bank firm survey 39 . The survey results are not public data. They are, however, available from the authors upon reasonable request and with permission of the World Bank. The figure data are available at https://figshare.com/projects/ Criticality_analysis_of_a_country_s_transport_network_via_an_agent-based_ supply-chain_model/91367. The simulation data generated by the model are available from the corresponding author upon reasonable requests.

Code availability
The code of the model is available at https://github.com/ccolon/ disrupt-supply-chain-model, a Git repository maintained by the corresponding author.


extfieldsofstudy:
Computer Science
Engineering

text:
Collaboratively Optimizing Power Scheduling and Mitigating Congestion using Local Pricing in a Receding Horizon Market

A distributed, hierarchical, market based approach is introduced to solve the economic dispatch problem. The approach requires only a minimal amount of information to be shared between a central market operator and the end-users. Price signals from the market operator are sent down to end-user device agents, which in turn respond with power schedules. Intermediate congestion agents make sure that local power constraints are satisfied and any potential congestion is avoided by adding local pricing differences. Our results show that in 20% of the evaluated scenarios the solutions are identical to the global optimum when perfect knowledge is available. In the other 80% the results are not significantly worse, while providing a higher level of scalability and increasing the consumer's privacy.

I. INTRODUCTION
Three trends set a challenge for future power grids. Firstly, the transition towards sustainable energy sources leads to more renewable energy, but also to a larger fraction of unpredictable and intermittent production. Secondly, the electrification of various systems such as transport (electric vehicles), heating (heat pumps), and in general an increase of electricityconsuming devices leads to a huge growth of power consumption. And thirdly, the distribution of energy generation leads to a very different pattern in the load of the power transmission grid, than it was designed for.
The control of a vast number of small power units, both consuming and producing, is extremely difficult to do completely top-down, so a centralized control strategy cannot be used [1]. At the same time the power infrastructure is aging, and was not built for the emerging pattern of distributed prosumers [2]. This is why we need self-organizing control algorithms to schedule the use of electric devices while taking into account constraints of the power distribution infrastructure. This aids distribution system operators (DSO) and transport system operators (TSO) to maintain power balance and make sure there is no congestion, i.e. the grid capacity is not overloaded.

II. PROBLEM STATEMENT
The problem at hand is a variation of the economic power dispatch problem [3], [4], where the operation of a set of generators is optimized, such that power is provided to consumers in the most cost-effective manner. Traditionally this problem would include controllable generators (power generation plants), constraining transmission resources (transformers, stations, cables), and end-consumers having a static load. Currently, with distributed energy resources, and demand response-the possibility to control the load of consumers using flexibility of smart devices-the problem changes significantly.
A formal definition of the problem is as follows: With all parameters defined in Table I, the objective function as defined in (1) is the same as for the traditional economic dispatch, which can be summarized as: to find a set of electrical powers x i for every device i ∈ N in the grid, that minimizes the sum of all costs J i (x i ). In the original dispatch problem defined, N defines only the producers. However, in our problem formulation (1) N denotes the full set of devices in the grid, including producers and consumers. This means that contrarily to the traditional economic dispatch problem, not only the generators are controllable, but also the flexible loads of end consumers.
Note that the costs in (1a) do not necessarily refer to the financial cost of a dispatch. Rather, this generic model only optimizes some social welfare, which defines a desirable outcome for all participants involved. Depending on the situation at hand, one could minimize the amount of greenhouse gasses emitted, or maximize the amount of renewable energy used. However, for the rest of this paper the costs are defined as the energy losses of the devices. By minimizing the amount of energy losses we attempt to find a dispatch that is both economic and sustainable.
Note that powers are indicated from the grid perspective. That means positive powers indicate power going from the grid to the device, and negative powers are from the device back to the grid. * Z indicates a type of agent, which can either be MO for market operator, CO for congestion, LOAD for a static load agent, PV for an agent with photovoltaic solar panels, or ST for a storage agent.
Let us discuss the constraints of (1), i.e. (1b-1e). The constraint (1b) states that the sum of all powers in the system has to meet a specific target Θ. In a balanced (island) grid this target has to be zero for every time slot, but in a connected grid this target is the contracted load with the transmission operator. Put differently, Θ is the net power input of the grid to the rest of the world. Constraint (1c) defines the limitations of the physical power grid; a topology matrix A specifies which device is connected to which grid components, and b is the rating of those components-the maximum amount of power that it can safely transmit.
Constraints (1d) and (1e) represent the constraints of the devices in the grid. Specifically, constraint (1d) states that a device cannot produce or consume more power than it physically can, which is represented by lower and upper bound, x min i and x max i , respectively. However, (1e) also takes into account another private constraint C. We refer to this constraint as private, as it only concerns the state of a single device, and its state may concern information that we should not require to share with other parties. E.g. an end-user should not need to share its intent to run his or her washing machine with the rest of the neighborhood. This constraint differs per device, and specifies device limitations such as a battery that cannot hold more than a certain amount of energy, and cannot discharge when already empty. Flexible loads also may have constraints considering the time at which they can turn on or off based on the consumer's settings. We will elaborate on these constraints, as well as the cost functions of the device agents in Section III-B.

A. Related Work
In existing studies, different strategies are used for energy management. Four different main categories are defined in [1] based on whether there is a distributed aspect of decisionmaking and whether there is one-or two-way communication.
One of the strategies defined in [1] uses Transactive Control, where distributed systems decide locally on their device management, using two-way communication in a marketbased control scheme. The authors compare this approach with traditional top-down switching, price-reaction and centralized optimization strategies and show that transactive control is capable of using the full flexibility potential of the smart grid devices, while maintaining the end-user privacy. This claim is consistent with earlier studies [5], [6] that have shown that using a market-based control in a multi-agent system can provide equally optimal results as a centrally optimized system, under certain conditions.
Multi-agent based methods are already getting attention in the smart grid domain due to many desirable properties: robustness, user-friendliness, attack resistance and scalability [7], [8]. The economic dispatch problem is very well suited to be represented using multi-agent systems [9]. There are many studies that use a multi-agent based approach to model and solve the problem, such as the two-way message passing agents using consensus algorithms to find an allocation that is optimal [10]. Other methods use a completely decentralized method involving reinforcement [11], utility maximization [12] or model predictive control [13].
In [14], a strategy is proposed to schedule the power consumption and production of generators and loads based on a method called Negotiated Predictive Dispatch. In this approach wind and conventional generators, as well as static and flexible loads, are controlled on the transmission level. Agents propose power schedules, which are aggregated by a central market operator, which then updates a price program using gradient descent, meaning the price of congested timeslots increase, and that of underused time slots decrease. In doing so the authors show that they are able to balance power production and consumption, while satisfying power grid constraints. However, this approach focuses on the transmission level, whereas we consider the distribution level power grid to be at a much more imminent risk of congestion. At the transmission level, the scale is much larger than at the distribution level, both geographically and considering the power levels involved. Moreover, the power grid has a very different topology at the transmission level. A major drawback of this approach when applied at the distribution level, is that using a global power price for congestion mitigation is not only unfair for agents in non-congested areas, but it would also be converging to a solution much too slowly. A similar approach to the negotiated dispatch is provided by [15] in which a gradient descent on the pricing is used, followed by a local optimization of agents.
Another approach is given by [16], in which power programs are proposed by device agents that are able to satisfy cluster constraints to reach a target energy consumption, while also staying within the limits of the grid. In their approach, called Profile Steering, agents are only motivated to accomplish the cluster goal, which is to consume or produce energy at a specific target amount, and will propose alternative power programs whenever constraints are violated. Proposals that reduce the constraint violations the most, are then selected as the new candidate programs. In our view, the problem of this approach lies in the lack of motivation for the participant to sacrifice private rewards for running an alternative program. The only objective is the cluster goal, which means that the price of having limited resources is paid by a select set of individuals that offer the most flexibility, or conversely, the profit is reaped by those who are able to maximally make use of remaining capacity.
The authors of [17] propose an approach based on DCOPs (Distributed Constraint Optimization Problem) to solve the economic dispatch problem, as well as the real-time demandresponse balancing problem. Their algorithm is able to find optimal solutions for a system of controllable generators and (predicted) loads, taking into account transmission network constraints. They show that finding the optimal solution using Dynamic DCOPs is possible, but their solution does not scale very well. A relaxed version of the problem, in which soft constraints may be (temporarily) violated, scales better, but still for relatively small time horizons, even considering their implemented solution on a GPU.
We note that related work on the topic of this chapter listed here is not exhaustive. For further detailed background, a comprehensive overview of different mechanisms for solving optimization problems in smart grids, particularly where demand response is involved, we refer to [18], [19], [20].

III. SELF-ORGANIZING ECONOMIC DISPATCH
Our contribution is a self-organizing method for coordinating the power scheduling of smart grid devices on the distribution level energy grid. This is a multi-agent based heuristic approach for find a solution to the economic dispatch problem as defined in (1). This multi-agent based approach allows for great scalability, while ensuring privacy and final control of the end-user.
At the distribution grid level, devices are typically consumer devices, such as PV-panels, household batteries, heat pumps, ventilation or air-conditioning units. The controllability, or flexibility of such devices is often limited, and bound by the device limitations and the user preferences. With increasing numbers of such devices, using a strictly "top-down" approach is intractable, since the search space grows exponentially with each added device. For this reason we use a hierarchical approach, in which the market operator delegates its task to intermediate "congestion" agents which then independently solve subproblems-which is to make sure that the total power throughput at their allocated point does not exceed a certain threshold. Making such subdivisions is justified, since in the power grid, transformers are effectively branches of the topology, and two nodes under different transformers are independent of one another; hence, transformers are the logical point to place congestion agents in the control hierarchy.
Our approach, denoted as Local Pricing Receding Horizon (LP-RH) is based on the economic incentive that end-user devices should have, to provide its owner with a service, in the most affordable way. Put differently, the system will use price differences to stimulate agents to schedule their power consumption or production in a balanced way, such that any grid constraints are satisfied. The overall scheme is simply to gather expected power programs from the connected devices, and then iteratively adjust energy prices to steer the agents into a certain power program. The method is explained in more detail in the following section, and is similar to the Negotiated Price Dispatch proposed by [14].

A. Local Pricing Receding Horizon
In order to create a power planning taking into account forecasts and/or predicted power programs of clients, we use a Receding Horizon (RH) approach, depicted in Fig. 1 which means that at any point we only take a fixed horizon of T Fig. 1. The principle of a receding horizon market is that every time step the time horizon for optimization shifts by one. Planned power programs (white) are turned into fixed contracts for the current time step (green), which determines the outcome of the algorithm. This figure is redrawn from [14].
repeat 5: for all j ∈ M i do 6: end for 8: 11: until ǫ ≤ ǫ max 12: end if 13: return x i program time units (PTU) into account. A PTU is typically 15 minutes or one hour, throughout this paper we will use PTU length of one hour so T = 24. However, none of the mentioned approaches are limited to this convention. When the algorithm has converged and a power program is found for the next T PTUs that satisfies all constraints, the first PTU becomes the "current" situation, and the projected power program becomes a contract. The time horizon now shifts by one, and the entire system starts again.
The Local Pricing Receding Horizon algorithm is shown as pseudocode in Algorithm 1; an example graph on which the algorithm could run is shown in Fig. 2. The algorithm describes how any agent finds a new power program x i to satisfy its local constraint C i . Every device agent in the grid does this by solving a local optimization problem (line 2), taking into account local constraints as explained in Section III-B.
Prices ρ are updated by the market operator agent and Simple example topology of a graph which could be running the LP-RH algorithm. One market operator (M ) is connected to two congestion agents (C i ), which are connected to a total of five device agents (D i ). The device agents can represent PV panels, consumer loads, storage agents, or any other leaf nodes in the power grid. Topologies where congestion agents are connected to more congestion agents are also possible.
the congestion agents in order to get a power program x satisfying the grid constraints. When the market operator or any congestion agent runs the LP-RH algorithm, all j ∈ M i , where M i are the immediate children of agent i, are requested to propose a power program based on the price ρ (line 6). Here the function is called recursively until all agents have determined a new power program. Note, that if any of the children are congestion agents, they determine their power program by forwarding the prices to their children and returning the sum of all received programs. Many agents will not know exactly how they will function in the future, as contexts might change, a user might behave differently than anticipated, or weather conditions may act up; this is why receding horizon method updates the time iteratively. With the sum of all power programs, the market operator can determine the error ǫ, which is the sum of all local constraint violations C for every t ∈ T . In line 9 the constraint C i is used to compute the local constraint violation and stored as an error variable ǫ. The value of ǫ is taken into account to adjust the prices in line 10. In the "adjustPrices" function, a gradient descent approach is used, which linearly interpolates the last two errors as a function of the price. We then choose the price at which the error is projected to reach zero. Note that this means we assume the reaction of the devices linearly depends on the prices, which will only hold under very specific circumstances. To overcome this issue, we iteratively repeat the process until ǫ ≤ ǫ max and assume that a non-linear response can be described as a series of linear pieces. The parameter ǫ max then represents an upper bound on the error, which we can use as a convergence criterion.

B. Agent Behavior
Agents in the system are characterized by the device that they have to assign a power program for. In the distributed case, every agent locally optimizes a local cost function J i , which is part of the global optimization problem (1). Also, local constraints C i have to be taken into account which are represented by constraint (1e). The behavior of the agents can be defined by three functions for the costs, the local constraint and the power program; an overview of this is shown in Table II. In our model we consider the following types of agents: 1) Market Operator: This is the root node of the tree (M in Fig. 4), as far as the local distribution grid concerns. In the physical grid, it corresponds to the transformer that connects the local low voltage (LV) grid to the medium voltage (MV) grid. We assume that there is no energy loss at the market operator, so J MO = 0.
Its goal is to find a solution to (1), and its private cost would be the same as the target constraint in (1b). The market operator runs Algorithm 1, using a deviation of the target profile Θ as an error, which is minimized by the algorithm. Hence, its local constraint is defined as wherex denotes the power program of the market operator.
Since the market operator is not a device in the grid itself, its power is the sum of the powers of its children, which in case of the market operator are all nodes in the cluster: The value of (3) can either be negative or positive, which respectively means either the total power production or the total consumption is too high. The market operator sets an initial price of ρ = 0.5 for all timeslots, and then uses Algorithm 1 to minimize the constraint value until it reaches zero, in order to satisfy the global constraint (1b).
2) Congestion Agent: This is an intermediate node on the grid tree (C in Fig. 4) connected to a parent node who is either the market operator or another congestion agent. It corresponds to a component in the grid where congestion might potentially occur, such as a transformer. Equivalently to the market operator, we assume that no energy is lost here, so: This agent has a constraint that aims to limit the power usage of that part of the grid, this corresponds to the constraint (1c). The congestion agent also uses Algorithm 1 to minimize the error of its local constraint, which is defined as Here β i is the power limit of the agent i, which means it is the maximum power throughput of the grid at the point i represents. Againx i denotes the power profile of the congestion agent, but is now equal to the sum of all devices under the current node: Similar to the constraint of the market operator in (3) the constraint value can become negative or positive, and is used to compute the error ǫ in Algorithm 1.

3) Load Agent:
This is an agent responsible for an uncontrollable load in the grid (represented by any device agent D in Fig. 4). This could be a consumer household, an office, street lighting, or any other non-flexible load, and thus only participates in the problem as part of constraint (1b). In the distributed system however, it is responsible for making a forecast of the power usage, and this forecast will be updated with more accurate information as the time horizon shifts.
The load agent has no attached cost in the global problem, and no need to locally compute any optimal behavior. This is equivalent to stating that its cost and constraint correspond to Its corresponding load profile x i is fixed to some profile that constrains the global problem (1). In our experiments its values are taken from real households as described in Section IV-B.

4) PV Agent:
The PV agent (any D in Fig. 4) has some flexibility to offer to the optimization function by allowing curtailment in reference to the expected generation. We assume that curtailment is binary, in that either the PV generates power as normal, or it is switched off and produces no power at all. Curtailing means that there is potential energy lost, and hence the cost function of a PV agent is defined as where x ′ i indicates the expected power, when not curtailing. This expected power is taken from the scenario, which will be detailed in the Section IV.
When reacting to prices in the distributed system, a decision is made in order to decide whether to curtail based on the price profile. If the operational running costs γ i of the PV is more than the power that would be generated by it, there is no point in running the generator (from an economic point of view). Hence, the local constraint and its corresponding decision rule of a PV agent can be written as For the PV agent it also holds that the corresponding expected power profile x ′ i determines the global problem (1). Its values in our experiments are taken from real PV panels as described in Section IV-B. 5) Storage Agent: A storage agent (again a leaf node D in Fig. 4) provides flexibility by allowing to store some energy in a local storage like a battery or a heat buffer. There are limits to the amount of energy that can be stored, either because of the physical limitations of the storage device, or because of the end-user settings. Moreover, a storage agent has some efficiency, which defines energy loss when energy is put into it, or out from it. This means that we can define the cost function of the storage agent as the energy lost during charging or discharging where This difference makes sure that the loss is correlated to the internal power of the battery when charging or discharging. I.e. x i defines the power at the grid side of the storage, when charging a lower power effectively charges the battery, and conversely when discharging a higher power is required to provide some power level to the grid. In order to define the constraints of the storage agent we must define the update function for the amount of energy e i stored as the cumulative sum of the powers Here λ i represents the leakage or self-discharge rate of the storage agent i, and e i0 is the stored energy at the start of the experiment. Let us now define the following constraints on the storage agent The minimum and maximum energy values are defined by e min i and e max i , respectively. For a storage device, the power limits denote the maximum charge and discharge rates. They are defined by x max i and x min i , respectively, in (1d). A special case of a storage device is a heat pump, which (electrically powered) heats a house in order to keep the temperature within comfortable levels. The heat pump only allows charging and only discharges through leakage, hence for a heat pump x min i = 0 and λ i > 0. When a storage agent has to update its expected power profile in line 2 of Algorithm 1, some response function (f (ρ) in Table II) is required, which is "economically sane" and has the following characteristics: • return x max i for low prices and x min i for high prices, • be a monotonically decreasing for increasing prices, • have a "plateau" of zero power response for some intermediate price (ρ = 0.5), which is wider for less efficient devices. The final characteristic allows more efficient devices to respond to subtle price change, and have less efficient devices respond to more extreme prices. This way devices with a higher efficiency are used first when flexibility is needed, and (when parameterized correctly) less efficient devices will only be used when required. In our implementation we chose the simplest response, where the agent linearly decreases its power from x max i to zero at ρ = η i /2. It then stays at zero symmetrically around ρ = 0.5 until ρ = 0.5/η i , and then decreases its power response linearly to x min i . This power/price relation is depicted in Fig. 3. An alternative strategy for the storage agent would be to make use out of any differences in the price, charging and discharging as soon as the price differences are large enough to overcome its efficiency. From a strictly economic perspective, this is the optimal strategy to maximize its own benefit. However, this leads to very "binary" behavior with minimal and maximal charging rates [12] and thus, little room for optimizing from the market operator and congestion agent. Therefore, a linear strategy is implemented as depicted in Fig. 3, allowing to solve the overall optimization problem (1).

IV. EXPERIMENT SETUP
The LP-RH algorithm was empirically evaluated by running a simulation of an LV grid with a set of realistic household load profiles and PV production profiles for a series of 24 PTUs.

A. Distribution Network Topology
For the topology of the network we use the European Low Voltage Test Feeder [21] network. This dataset is used to benchmark power and energy algorithms on realistic European distribution networks. In this paper we superimposed six points on the topology, where we monitor and mitigate any potential congestion. These points are strategically chosen to separate the problem into independent subproblems. The resulting topology with the congestion points are shown in Fig. 4.

B. Household Load and PV Profiles
The household consumption and production profiles are taken from a pilot study [22], in which 92 residential consumers were monitored over the course of a year (from March through November 2018). The data was preprocessed such that we have separated information on the consumption of houses, and of the PV installations. Data is anonymized and randomized per month, so that we can select data from any specific month for a base load of a household, or a residential PV installation.
In the experiment the 54 households from Fig. 4 were assigned a random instance of the load and PV profiles from the same month (i.e. all households were equipped with PV panels). The daily load consumption varied between 4.98 kWh and 29.39 kWh, and the total PV production varied between 826 Wh and 18.8 kWh. Furthermore, 16 randomly selected Market operator Congestion agent Household Fig. 4. The topology of the IEEE LV feeder network used for the simulation. The figure depicts the connections between household consumers, the intermediate congestion agents, and the root market operator. The "filled" markers represent the agents from Section VI. households were assigned a household battery, and again 16 were chosen to have a heat pump installation.
Every household and PV installation in the simulation would select a random profile from the dataset, which was used as its power profile. The objective Θ was set to a total net consumption of the LV grid using the mean total power consumption of the houses including PV production. Choosing the target profile Θ in this way corresponds to a situation in which the energy provider of the simulated neighborhood would agree to a contract for the average behavior of the households, and subsequently attempts to use flexibility to account for any deviations from the normal.
The batteries were all dimensioned with storage capacities of e max = 10.8 kWh, and maximum charge and discharge rates of x max = 4000 W and x min = −4000 W, respectively. The batteries charging efficiencies were all set at η = 0.9. The heat pumps were estimated to have a working energy capacity of e max = 2 kWh, this means the difference between the thermal energy capacity of the house at the minimum and maximum comfortable user temperature is 2ψ kWh, where ψ = 4.0 is the coefficient of performance of the heat pump. Then, the heat pumps have a x max = 1600 W and x min = 0 W, an efficiency of η = 1 and a constant leakage rate of λ = 360 W. Finally, to ensure convergence, the maximal error value is set to ǫ max = 10 −3 in the experiments for this paper

C. Forecast Uncertainty
The predicted load and production power profilesx of the load and PV agents were generated by taking average profiles of the complete dataset. These average profiles are considered as taken from historic data and hence, provide a ground for predicting the power program of future PTUs. When agent i determines its power prediction x ′ i , it will compute a weighted average between its assigned power profile x i and the average profilex, such that: such that at t = 1 the prediction equals the selected profile x ′ it . At t = T , the prediction is simply the average powerx t .

V. CENTRALIZED SOLVER
In order to address the performance of the LP-RH algorithm, a centralized optimization approach is also introduced to provide lower bounds to its results. A mixed integer linear program (MILP) solver was used to find these bounds for the problem stated in (1). The centralized optimization approach considers the exact same scenario that was solved by the decentralized algorithm; the energy loss is minimized and the same set of constraints apply. A fundamental difference lies in the availability of information. Whereas detailed information is only shared locally in the LP-RH algorithm, the centralized optimization approach assumes complete knowledge of the current system state for the decision-maker; i.e., no limits are imposed on the spacial flow of information within the network. With these features in mind, two versions of the MILP were formulated: Receding Horizon Centralized Solver and Perfect Information Centralized Solver.

A. Receding Horizon Centralized Solver
The receding horizon centralized solver (RHCS) is the most similar to the LP-RH algorithm. It uses a receding horizon approach (as depicted in Fig. 1) to find an ideal dispatch solving the consecutive sub-problems. Uncertainty about future device states is again simulated by (17). Because this solver has perfect information at the current state for each iteration of the receding horizon, its solution represents a lower bound on the solution found by the LP-RH algorithm.

B. Perfect Information Centralized Solver
In the perfect information centralized solver (PICS), a single centralized optimization problem is solved for the complete dispatch. In addition to perfect spacial information, this version of the MILP also has perfect temporal information about the complete system state; this means that no uncertainty about future states is simulated. Referring back to (17), this is equivalent to setting α = 1, which results in perfect predictions for each profile. The solution of this MILP represents an absolute lower bound for problem (1).

VI. RESULTS
In this section the results of the described experiments are shown. Before comparing the performance of the different solvers, which is done in Section VI-A, the behavior of the LP-RH algorithm is shown in this section. The graphs in this section show examples of single simulation runs, demonstrating the behavior of the different agent types. Powers are shown as power consumption, this means a net consumption is shown as positive power, and conversely negative powers indicates a net power production. In Fig. 5 the final result of the market operator is shown, where the market operator has found a price profile such that the target profile is met exactly. Fig. 6 shows the power profile of one of the three congestion agents that is directly connected to the market operator. Its power congestion limits are set such that in the peak moments of the day there is some congestion expected. This results in a price difference shown in the power program, as around the peak PV production (t = 12, 13, 14) the local price is slightly lower than the market price, leading to a lower net power production. Similarly, at the end of the day (t = 19, 20, 21, 22) a positive power congestion was mitigated by increasing the price.
The power program of a heat pump agent is shown in Fig. 7. This heat pump is connected directly to the congestion agent from Fig. 6, hence its price profile should be identical. What is most obvious in this graph is that the high price at t = 12 leads to a zero power consumption, and quickly after that, the power consumption rises in order to maintain comfortable temperatures. Again, around t = 19, 20, 21 the relatively high prices lead to a power consumption of zero, which was apparently feasible because of the high power consumption leading up to it.
Finally, Fig. 8 shows the total power programs of all devices of the four classes in this experiment summed up. In this figure, the power profiles of the loads and the PVs are the direct result of the chosen profiles, and are the input for problem (1). We can see that for the majority of the experiment, all used flexibility is from the heat pumps with the high efficiency. Only at times of the congestion will the less efficient battery agents be used.

A. Comparison with Central Solvers
In a validation experiment 105 random problem instances (15 permutations for seven months) were created, and solved by the three different algorithms. In Fig. 9 the results are shown for the feasible instances. A solution is considered  feasible if all three solvers were able to find a solution that satisfies all constraints. The LP-RH algorithm did not find a correct solution for 26 problem instances, 11 of which were found to be overconstrained according to RHCS. For the feasible instances the LP-RH algorithm found solutions that were not significantly worse than the PICS, and in 21 instances found a solution with the exact same cost. In 27 instances LP-RH found a solution that was equally well as the RHCS or even slightly better-this seems to contradict the initial statement that RHCS acts as a lower bound for LP-RH; however, this is due to the way both solvers deal with uncertainty. The LP-RH algorithm allows the congestion agents to violate power constraints for future PTUs, as long as eventually they are not actually congested due to imperfect forecasts according to (17). The RHCS algorithm does not allow power constraints to be violated at any point in the future, Power loss (kW) Fig. 9. Compared with centralized optimization solvers, the LP-RH algorithm performs equally well. In this box plot the median is shown as a line, the boxes indicate the 25% and 75% percentiles, and the tails are capped at a maximum length of the box width; outliers are drawn separately. and hence might react too strict when a future congestion is mistakenly predicted. In the set of infeasible problems, i.e. problems that do not have a valid solution satisfying all constraints, LP-RH does find solutions that minimize the power loss. However, since these solutions either temporarily overload congestion agents (where x a gen > β i ), or does not match the target profile Θ exactly, they are not a fair comparison, since they do not strictly solve (1). In a separate run the problems were relaxed, by increasing the power rating of the congestion agents. This resulted in the LP-RH not being able to find a feasible solution in only 14 problem instances, but the results were otherwise very similar to the ones reported here in Fig. 9.

VII. CONCLUSIONS
We have introduced an algorithm for self-organizing smart grids, by solving the economic dispatch problem using a decentralized market based approach with local pricing. The LP-RH algorithm using a hierarchical approach was shown to be able to solve the problem using a fairly simple interaction scheme in which pricing information is sent down the hierarchy tree, and planned or forecasted power programs are sent back up. Using a gradient descent approach, the market operator is capable of tuning the pricing to find feasible solutions to minimize the power losses in the grid.
In our experiments we found that in 20% of the problems LP-RH did not perform any worse than a perfect-information centralized solver. In the other 80% our algorithm did not perform significantly worse. The benefit of LP-RH over a centralized solver are in the robustness and scalability of the solution, as well as in the preserved privacy of the endconsumers.
In the implementation of the response of the storage agent, we intentionally did not choose to respond with an optimal power program given the price signal. Particularly, when a high price is expected in the future, the agent will not "proactively" charge to avoid having to charge later, or vice versa. This behavior could be implemented at the agent quite easily using a dynamic programming approach, but it would lead to very extreme behavior, e.g. very binary behavior of charging or not-charging at full capacity even for small price differences. This binary behavior is hard to deal with in the rest of the hierarchical tree, and does not lead to any problems per se, but might be improved upon in a future continuation of this work.
Other variations of the problem may include other device types, for instance time-shiftable devices such as washing machines or dishwashers. Also, using electric vehicles (EV) as an additional type of agent, providing energy flexibility is a very interesting extension, which will undoubtedly lead to other complications because of their high power ratings. Finally, an integration with a real time balancing algorithm such as [23] would be very fruitful to complete the needs of the future smart grid.


Urban Transportation: Performance and Problems (Case Study: Route of ABG, CKL, and AT)

ABG, CKL, and AT are public transportation (called angkot) with the longest route in Malang. The length of the routes of ABG, CKL, and AT is 26m, 22Km, and 18Km. This research aims to determine the performance of three public transportation and to provide a new route recommendation. The operational performance is assessed from headway, load factor, and travel time. The service performance is assessed from safety, comfort, affordability, and regularity. Result of the operational performance shows that the load factor of three public transportation has met the specified standard, the headway and the travel time of ABG and AT has met the standards, however, CKL’s headway and travel time have not met the standard. The ImportancePerformance Analysis (IPA) shows that there are some service attributes which need to be improved. The recommendations that can be given are shortened the trajectory of ABG and extend the track for CKL and AT.

I. INTRODUCTION
Transportation is a type of service that can be used to move goods from an area (origin) to another area (destination) [1], [2]. Transportation requests will arise if there is a purpose in the request, for example, the desire for recreation, the desire to work or go to school, the desire to shop and so on [3]. Transportation requests are inseparable from the needs and availability of transportation that serves the movement of people and goods.
One of the urban transportations that can be used by society is city transportation [4]. City transportation (angkot) is a facility that serves the movement of both passengers and goods in an urban area. Angkot is in the form of motorized vehicles and used by the community with certain fares collected directly or indirectly [5]. Angkot has a predetermined route and different travel times in serving community movements.
Malang has 25 public transportation routes with the number of vehicles in operation is 1153 vehicles [6]. Each angkot has different routes and destinations. The optimal city transportation route can be seen from the distance and travel time of a trip by paying attention to the density of the road at a certain time [7], [8]. The density of the road at a certain time can affect the travel time of an angkot from its origin to the destination. If the traffic volume exceeds the existing road capacity, it will cause transportation problems which are called congestion [9]. ABG, CKL, and AT are angkot with the longest routes compared to other angkot in Malang City. The length of the ABG route is 26 Km which stretches from the northern part of Malang City to the south. The length of the CKL route is 22 km that stretches from the east to the west of Malang City. Whereas AT has a route length of 18 Km that stretches from the north to the west of Malang City. The three angkot pass through the road sections which often experience congestion at certain times (peak time), namely Sunandar Priyo Sudarmo Street, Kol. Sugiono Street, Gatot Subroto Street, Ranu Grati Street, and Panji Suroso Street [6]. Those roads are not only passed by ABG, CKL, and AT but also by many other angkot with different routes so that overlapping routes occur. The overlapping route can affect service performance and operational performance of ABG, CKL, and AT.
The operational performance of ABG, CKL, and AT have not met the standards set by the Director-General of Land Transportation in 2002, i.e. the standard for load factor is 70% during peak hours and non-busy hours and the headway set is 2-5 minutes during peak hours and 5-10 minutes when the hours are not busy. The load factor for ABG routes is only 31% with 19 minutes of headway [10]. The CKL route load factor is 34% with 18 minutes headway while the load factor AT route is 43% with 27 minutes headway [6].
The objectives of this study are twofold. The first is to determine the performance of public transportation in Malang City, especially the ABG, CKL, and AT routes, by assessing the operational performance of public transportation and service performance. The second is to provide recommendations for optimal city transport route in Malang City based on the trip generation of Malang residents.
II. RESEARCH METHODOLOGY Stratified Ramdom Sampling is used to calculate the samples. The sample size is determined using the Slovin formula [11].
where, n = number of sample N = population e = significance level =10% Sample for each angkot can be calculated using (2). Summary of research sample is shown in TABLE I. n = (angkot/population) x Sample  Headway analysis is used to determine the difference in time needed between two city transports with the same route [4], [13].

Headway (H) = 60/Frequency 
Load Factor Analysis is used to compare the number of passengers in city transportation with the transport capacity [4], [12].
Load Factor (LF) = (No. of Passengers/Capacity)x100%  Travel Time Analysis aims to determine the effectiveness of public transport by calculating the travel time of public transport trips from origin to destination [4], [13].

W = T/J 
where, W = Travel time (minute/km) J = Distance between segment (km) T = Angkot travel time (minute) Importance-Performance Analysis (IPA) used to measure the level of importance along with the satisfaction of angkot customers towards the services of security, safety, comfort, affordability, and regularity of angkot in Malang City [14]. The level of conformity is calculated using: where, Tki = Level of conformity Xi = Performance score Yi = Importance score The average score for each variable was plotted in the two-dimensional state space to create the IPA diagram. The vertical axis illustrates the importance score, while the performance score is labeled by the horizontal axis. The IPA diagram consists of four quadrants, i.e., concentrate here, keep up with the good work, low priority, and possibly overkill. The first-quadrant, i.e., concentrate here, which is located in the north-west corner is the one with low performance but importantly perceived by the customers, therefore the company should invest more to improve these attributes so the customers will be delighted. The secondquadrant is kept up with the good work. It is the one that is considered as important and the customers are fond of the performance of the service. The third-quadrant is a low priority. The attribute belongs here are performing well but customers perceive them as less important when compared with other attributes. The last or the fourth-quadrant is considered less important by the customers and felt too excessive so that need to be reduced due to the excessive investment. Trip Generation Analysis aims to estimate the amount of movement from the origin zone to the destination zone [3], [7]. Malang City is divided into 27 zones spread in each sub-district.

A. Passenger Profile
The participants of this survey were required to have been experienced in doing transactions with the object of the study. The potential participants were first approached and asked if they agreed to participate in the survey. The profile of the passengers is shown in TABLE II.    Fig. 1. The item statements belong to the first quadrant are smoking prohibition (6), waiting time (9), and route conformity (13). It means that the service provider (ABG) needs to pay attention to those item statements and invest more to gain satisfactory from the customers. The statements belong to the second quadrant are exit and entrance (3), load factor (4), windows (5), fares (7), and travel time (11). Since the second quadrant shows the attributes that have high performance and high importance values, the ABG needs to keep up its good work to maintain a satisfactory level from the customers. The statements belong to the third quadrant are lighting (1), window film usage (2), and headway (10). The ABG angkot does not need to invest more on those item statements because the customers do not consider them as important. The last, the statements that belong to the fourth quadrant are angkot availability (8) and arrival schedule (12). These attributes need to be reduced due to excessive investment.
The IPA diagram of the result of the CKL angkot study is shown in Fig. 1. The item statements belong to the first quadrant are waiting time (9), headway (10), travel time (11), and arrival schedule (12). It means that the service provider (CKL angkot) needs to pay attention to those item statements and invest more to gain satisfactory from the customers. The statements belong to the second quadrant are exit and entrance (3), load factor (4), windows (5), fares (7), angkot availability (8), and route conformity (13). Since the second quadrant shows the attributes that have high performance and high importance values, the CKL angkot needs to keep up its good work to maintain the satisfactory level from the customers. The statements belong to the third quadrant are lighting (1), window film usage (2), and smoking prohibition (6). The CKL angkot does not need to invest more on those item statements because the customers do not consider them as important. The last, there are no statements that belong to the fourth quadrant.
The IPA diagram of the result of the AT angkot study is shown in Fig. 1. The item statements belong to the first quadrant are waiting time (9), headway (10), arrival schedule (12), and route conformity (13). It means that the service provider (AT angkot) needs to pay attention to those item statements and invest more to gain satisfactory from the customers. The statements belong to the second quadrant are window film usage (2), load factor (4), fares (7), angkot availability (8), and travel time (11). Since the second quadrant shows the attributes that have high performance and high importance values, the CKL angkot needs to keep up its good work to maintain the satisfactory level from the customers. The statements belong to the third quadrant are lighting (1), exit and entrance (3), windows (5), and smoking prohibition (6). The AT angkot does not need to invest more on those item statements because the customers do not consider them as important. The last, there are no statements that belong to the fourth quadrant.

IV. CONCLUSIONS
From the analysis of operational performance using headway, load factor, and travel time analysis, it is known that the ABG and AT headways are in line with service standards while the CKL headway still has not met headway standards. The load factor for this angkot is still below 70%, it means that the load factors are not under the standard. While the travel time for these three angkot has met the standard except for the travel time for CKL in the direction of Landungsari -Cemoro Kandang Terminal on weekends and weekdays. From the analysis of service performance using IPA, the attributes of ABG angkot services that need to be improved are smoking prohibition, waiting time, and route conformity. The attributes of CKL public transportation services that need to be updated are waiting time, the time between, travel time, and arrival schedule. While the attributes of AT angkot services that need to be improved are waiting time, the time between, arrival schedule, and route conformity. The route recommendations for the ABG is to shorten the route from a length of 26 Km to 22 Km. Where as the other two angkot routes are extended from 22 Km to 22.3 Km for CKL angkot and from 18 Km to 18.5 Km for AT angkot.


Entrepreneurship capability in industrial engineering and its effect on competitiveness

Entrepreneurship and competitiveness are interesting topics in regional economic development. The study aims to examine the effect of entrepreneurship on competitiveness in the leather industry. The study was applied to the leather tanning industry in Indonesia. The research method is based on a survey of 42 leather tanning companies, data analysis using the Partial Least Square. The research findings the effect of Entrepreneurship on Competitiveness has a low correlation. The cause of the low level of Entrepreneurship in SMIs is mainly in achievement motivation, leadership and business orientation factors. Therefore it is suggested to employers and the government to encourage entrepreneurial skills through entrepreneurship education, facilitation of capital and marketing and building partnerships in the leather industry.

Entrepreneurship does not only apply to large companies, but also to small and medium enterprise. In this case the entrepreneurial relationship and competitiveness are interesting to study, especially in the SMIs manufacturing industry in Indonesia. Therefore, the research aims to measure the effect of entrepreneurship on competitiveness in Small Medium Industries (SMIs) so that this research is useful for entrepreneurs and the Government so that the performance and competitiveness of the industry can be improved.

Literature review
Porter discusses the relationship of a country's welfare reflected in industrial competitiveness, the relationship of state competitiveness with industrial competitiveness is illustrated in the Diamond model [4]. Economic competitiveness in an area is influenced by the structure of competition between companies in an industrial cluster [9]. While other researchers call competitiveness as success in achieving goals, where success is determined by the capacity of companies in the industry [5]. But researchers have not agreed on the measurement of competitiveness because each level of the organization has dimensions with different indicators [10].
There are differences in the measurement of competitiveness in the manufacturing industry. The characteristics of the manufacturing industry have raw material entities and machine facilities with clear stages and process flows [11]. This is as a manufacturing system including the flow of raw materials, the production process so that it becomes a product to be submitted to consumers [12]. Competitiveness is influenced by the interaction of internal and external factors as competitiveness driven. Internal driven factors are specified as Managerial and Capability of technology, while external driven factors are the Industrial Structure and the Role of Government [13].
Referring to the various opinions above, understanding Competitiveness is the ability to achieve future goals successfully and superior to competitors. Furthermore, in the manufacturing industry, competitiveness consists of three dimensions, namely Resources, Operational Processes and Performance [14]. Dimension grouping refers to the stages of the process in a manufacturing system, namely: resources (input), processing (process) and product output (output).
Moderate understanding of entrepreneurship is one's ability to manage, manage and take business risks [15]. While in Indonesian the term entrepreneur, consisting of the words 'wira' which means brave, mighty, and 'business' which means doing business, so entrepreneurship means someone who is brave, able to run a business to achieve success independently. The term entrepreneurial or entrepreneurship is used interchangeably as an activity in running a business, namely as an entrepreneur.
In economic dynamics challenges and opportunities and risks arise in the midst of uncertainty in business. Entrepreneurship related activities to find, evaluate and develop opportunities, business plans, resource allocation and business control [8]. The series of company activities is controlled by the leadership as an entrepreneur. Hunter & Lean mentioned the role of company leaders in business capability is reflected through the character of Entrepreneurial Leadership [16].
Matsuno et al. referred to entrepreneurship as determination, enthusiasm to improve the company's business performance by taking into account the orientation of meeting market needs [17]. Jones and George added that entrepreneurial abilities include experience, risk taking, high self-esteem, a sense of ability, the need for achievement [18]. While Lee & Hsieh called the ability of entrepreneurship as the ability of marketing, innovation, to achieve sustainable excellence [19].
Referring to the opinion above, entrepreneurial abilities include enthusiasm, courage to take risks and achievement motivation, manage resources with organizational leadership to gain profits and achieve excellence on an ongoing basis. In this study the entrepreneurial ability variable has three dimensions, namely: Entrepreneurial Spirit, Leadership and Goal Orientation.
The linkages of entrepreneurship to competitiveness were examined by Buckley et al. and Man et al. [20,21]. They analysed the effect of entrepreneurship on competitiveness in SMI companies. Entrepreneurial ability translates as leadership and managerial. Furthermore Man et al. examined entrepreneurial competencies in small and medium industries (SMIs) and socio-cultural influences on the formation of entrepreneurial competencies [22]. He found that entrepreneurial competence was a high characteristic related to company performance.
This study discusses the effect of entrepreneurship on competitiveness in the leather industry in Indonesia. The entrepreneurial ability of SMIs is also reviewed based on the conditions and scope of competition and organizational capability. The pattern of the relationship between Entrepreneurship Capability to Competitiveness was also examined by Bulis-Skapars [23]. He found that the effect of Entrepreneurship on industrial competitiveness would fail if it did not pay attention to environmental factors in the region. The study discusses the effect of entrepreneurship on competitiveness in the manufacturing industry, namely the tannery industry.

Research methodology
This Research aims to examine the effect of entrepreneurship on competitiveness is carried out in the tanning industry. The structural model with formative patterns and validated empirically in the manufacturing industry. The population consists of 300 consisting of 50 formal companies and 250 informal firms of the SMI's leather tanning industry in Garut, Indonesia. Data collecting on purposive sampling (non-probability sampling) to 50 formal companies, with the number of questionnaires collected 42 companies or a response rate of 84%. Respondents are entrepreneurs as owners or factory managers, while competitiveness is measured through self-assessment. The questionnaire contains the respondents' responses as a Likert scale which has the smallest range of values = 1 and the largest = 5. Data processing is displayed descriptively, and analysed data using Partial Least Square (PLS).

Result and discussion
The results of analysis and processing data of Entrepreneurship (EP) influence Industrial Competitiveness (IC), shown below:  Table 2 and Figure 1. Path analysis shows that the entrepreneurship variable has a coefficient ß1 = 0.674 with a determination (R2) of 0.454. Effect of Entrepreneurship on competitiveness (R2) 0.25-0.5 as a weak correlation. Furthermore, the model capability test (predictive relevance) with Stone-Geisser's obtained Q2 = 0.213, the criteria of Q2> 0.15 indicate the model's capability as moderate  Figure 1 shows the structural model reflects relationship between variables. The coefficient of determination (R 2 ) = 0.454 this means that the three variables can explain competitiveness variance of 45.4%. The test results also showed that the entrepreneurship variable was strongly reflected by the Leadership dimension of 0.973, while the Entrepreneurship (EP) dimension was 0.526 and the weakest was the Goal Orientation dimension (GO) of only 0.247. The model capability test with the Stone-Geissler's parameter denoted by Q 2 . Capability Test reflects the ability of the model to predict relationships between variables (predictive relevance). PLS analysis with blindfolding procedure shows the value of Q 2 = 0.276, the category of value of Q 2 > 0.15 shows that the model has great capabilities so that the model is stated as moderate.   Table 2 shows that variable has a positive with coefficient of ß1 ≠ 0, t-value > 1.96 so H0 is rejected. Thus the alternative hypothesis (H1) is accepted and stated that Entrepreneurship have a positive effect on the competitiveness. The results of analysis show the overall model is fit but the correlation is weak.
Related to the theory, why is the model entrepreneurship has a low correlation with competitiveness? This explanation refers to a descriptive analysis of the educational background of employers in the leather industry. In theory, education has strong influence on entrepreneurship [24]. While, the education level of leather entrepreneurs is low. The proportion of entrepreneur who have graduated from college is very small, only 14,3%. The greater proportion each: have a high school education of 40.5%; 31% of junior high school and 9.5% graduated from elementary school, and the others 4.8%. The level of education of workers is generally a senior high school.

Conclusion
The results found the effect of entrepreneurship on competitiveness has a low correlation. The reason for the low influence of entrepreneurship on competitiveness is: there is a weakness in the dimensions of entrepreneurship. There are weaknesses in the organizational leadership dimension, the dimension of goal orientation, and the entrepreneur's educational background. The dimensions of organizational leadership with creative and innovative indicators showing low scores. The dimension of Goal Orientation with sales achievement indicators and profit achievement indicators also showing low scores. While overall the entrepreneurship is moderate, meaning the capability of entrepreneurship has not been categorized as good.
Based on their educational background, the majority of SMIs entrepreneurs have low levels of education. Therefore it is necessary to improve education so that they able to improve the entrepreneurial capability and competitiveness of SMIs in the future. Related to the limitations of the scope of research, in addition to considering the level of education, we suggest conducting research in other industrial sectors with a larger sample size. And need to include other variables such as regional economic conditions that tend to decline, Related to the limitations of the scope of research, in addition to considering the level of education, we recommend conducting research in other industrial sectors. As well as the need to include other variables such as regional economic conditions that tend to decline, in order to obtain comprehensive and useful research results in improving the performance and competitiveness of the industry in general