Argonne National Laboratory
Idaho National Laboratory
Oak Ridge National Laboratory
Pacific Northwest National Laboratory
Methods: Data-driven stochastic optimization, Optimization under endogenous uncertainty, Bi-level optimization, Mixed-integer programming, Decomposition algorithms, Statistical-learning-based optimization, Active learning, Reinforcement learning, Simulation, Machine learning, Parallel computing
Applications: Critical infrastructure protection, Sustainable and electrified transportation network design (e.g., aerial drones, electric vehicles), Security of cyber-physical systems (e.g., smart manufacturing, smart grids), Nuclear alloy design, Electricity markets and Power systems operations with renewable energy integration, Power grid resiliency, Vehicle routing, Cyber and information network security, Wireless network security, Disaster mitigation, and Supply chain resiliency and risk mitigation.
In many real-life decision-making problems, practitioners are often confronted with uncertainty in the underlying system parameters. Moreover, in solving problems under uncertainty, analysis of historical/survey data can inform the decision-makers about uncertainty such as decision-makers can estimate the probability distributions of the uncertain parameters analyzing data. Also, decision-making problems often involve multiple competing players, requiring the decision-makers to account for other competing players' objectives in best accomplishing their objectives. My academic interest is to develop mathematical optimization frameworks for problems having uncertainty and capturing the interaction of multiple competing players. Specifically, I develop data-driven stochastic programs and bi-level programs to model these problems as well as develop customized decomposition and approximation algorithms utilizing the problem structure to solve those optimization problems.
Besides optimization, practitioners often deal with problems that require analyzing big data to identify hidden patterns such as detecting malicious cyber nodes in a large computer network. My research interest also lies in developing/implementing supervised and unsupervised machine learning algorithms for big data analytics utilizing high-performance computing environments.
We can utilize machine learning models in optimizing problems where the traditional derivative-based optimization methods are not applicable. I also utilize statistical-learning-based optimization methods, such as Bayesian optimization, to solve decision-making problems where derivative information is not available or it is computationally very challenging to estimate.
My Ph.D. dissertation is focused on developing new data-driven stochastic programming models, mathematical reformulations to linearize the original nonconvex models, and designing customized decomposition and approximation algorithms for solving optimization problems under endogenous (decision-dependent) uncertainty. Unlike traditional stochastic programs, where uncertainty is usually modeled using scenarios governed by a known probability distribution, in some application areas, such as retrofitting systems under imperfect and multi-level effects of protection resources on component reliability, it is beneficial to model the probability distributions of the uncertain parameters as a function of the decision variables. The stochastic program modeling this decision–dependent uncertainty is known as stochastic programming with endogenous uncertainty (SPEU).
In a two-stage stochastic program (as shown in Figure 1), scenarios are generated to represent the realizations of uncertain parameters following a probability distribution, where each scenario represents a particular realization of the uncertain parameter. The optimization problem is decomposed into two-stage, where the first-stage decisions (x) are made without realizing the uncertainty. But the second-stage decisions (a.k.a., recourse decisions) are made after realizing the uncertainty, meaning the second-stage model is solved for each particular scenario after realizing the uncertainty. In optimization under endogenous uncertainty, decisions alter the probability distributions of uncertain parameters. For example, if a decision-maker invests more resources in a network component (e.g., better quality material) then the survival probability of that component will be higher when exposed to disruption, making the survival probability distribution to be a function of investment decision.
Therefore, as the probability distribution of uncertain parameters is a function of the decision variables in an SPEU program (Right side in Figure 1), the scenario probability is also a function of the decision variables (x). Therefore, in the first-stage objective function (Right side model in Figure 1), scenario probability is a function of the first-stage decision variables (x). But, this is not the case in the stochastic program with exogenous uncertainty (Left side model in Figure 1), where scenario probabilities are fixed, as the realizations (as well as the corresponding probability distributions) of uncertain parameters do not depend on the decisions (i.e., exogenous to the model). However, in this stochastic program with exogenous uncertainty, the first-stage decisions affect the recourse decisions by appearing in the second-stage model as a parameter (Left side model in Figure 1).
This SPEU framework has many application areas including but not limited to critical infrastructure protection, disaster preparedness, and mitigation, supply chain risk mitigation, power generation capacity expansion planning. Solving problems having a natural decision-dependent uncertainty structure using a traditional (exogenous) stochastic program could be computationally very difficult as the exogenous stochastic program would need a larger number of scenarios than the SPEU framework to capture the decision-dependent uncertainty property. Moreover, modeling and solving problems as SPEU can exploit the underlying decision-dependent uncertainty property of the problems and can yield better quality solutions and more insights into the problem. But, this SPEU modeling framework has not been utilized in many application areas to date as well as most of the studies on SPEU approached the problem with problem-specific heuristic methods and commercial nonlinear solvers. But, heuristic methods cannot provide optimal solutions, and commercial optimization solvers/software are not available to industry decision-makers always. Moreover, directly using optimization solvers/software is not computationally efficient for solving mathematical models with difficult structures (e.g., nonconvex stochastic mixed-integer nonlinear). Also, very few studies analyzed theoretical properties (e.g., convexity, submodularity, etc.) of the SPEU framework and limited them to only binomial probability distributions. But, there are other widely used probability distributions in many areas. It is possible to provide good quality approximate solutions to large-scale problems by taking advantage of some theoretical properties. Therefore, there is a need to make advances in the SPEU modeling and solution approaches. Until this need is met, decision-makers in various application areas cannot obtain good quality solutions and valuable insights into their decision-making problems, resulting in systems that underperform and have excessive cost.
Therefore, in my Ph.D. research, I have developed novel data-driven stochastic programming with endogenous uncertainty (SPEU) models for solving new problems having endogenous (decision-dependent) uncertainty that has not been studied to date. These developed mathematical models incorporate supervised machine learning models into the stochastic programming framework to estimate the probability distributions of uncertain parameters by analyzing historical/survey data. As in the SPEU models, scenario probabilities are a function of the decision variables of the optimization routine, the models are naturally nonconvex. This nonconvexity makes it difficult to solve these SPEU models compared to the traditional stochastic programming models. I developed new mathematical reformulations to linearize the developed nonconvex nonlinear models to reduce computational complexity. I designed customized decomposition algorithms to solve the reformulated models for realistic-sized problem instances to optimality. I also incorporated some new valid inequalities to improve the computational efficiency of the developed algorithm to outperform the existing solvers in computation time.
In my dissertation, I studied the submodularity property of the SPEU modeling framework and mathematically proved that the objective function of the stochastic program is submodular when the uncertain parameters are exponentially distributed. Taking advantage of this submodularity property, I implemented an efficient approximation algorithm capable of solving large-scale decision-making problems under endogenous uncertainty much faster than the exact algorithms with an approximation guarantee to the solution quality. Unlike the heuristic algorithms, approximation algorithms provide a worst-case performance guarantee to the solutions. In my Ph.D. research, I demonstrated the benefits of the developed methodologies in the application areas of infrastructure protection and reliable transportation network design, disaster mitigation, reliability optimization of complex systems, and supply chain risk mitigation.
Computer and Information Network Security
My foray into cyber-security research was through my M.S. research at Mississippi State University, (1) cyber network security problems from the perspective of a cyber network owner (defender) who seeks to efficiently utilize the limited defense resources to deploy security countermeasures (e.g., firewalls, intrusion detection/prevention systems) to protect critical assets from cyber-attacks. In reality, cyber systems are vulnerable to multiple types of adversaries (attackers), each with different capabilities (e.g., skills, resources, etc.). Additionally, network owner faces uncertainty about the attackers’ capabilities, resources, and probabilities of attack success while making strategic decisions against them. In the presence of these uncertainties about attackers, it is crucial to find an optimal security countermeasures deployment decision that performs well under any particular attacker’s attributes (e.g., capabilities, resources, attack success probabilities).
To provide decision support to the cyber network owners, I used the concept of an attack graph that maps an organization’s cyber network topology to potential attack paths through which attackers can breach the critical assets of the organization. A node in the attack graph represents a security condition or an asset, whereas an arc represents an attacker’s action (attack). The goal of the network owners is to deploy countermeasures in the network that interdict (remove) the corresponding attack paths from the network, meaning that the attacker cannot use that path. Figure 5 demonstrates how attack graphs can be used to represent the potential attack paths for an example cyber system.
Figure 5: Attack graph representation of an example cyber system
In one of the projects, I developed a two-stage stochastic mixed-integer programming model based on the attack graph concept to minimize the expected loss from cyber attacks by optimally deploying security countermeasures under uncertainty in attackers’ attack success probabilities [1]. We implemented a Sample Average Approximation (SAA) algorithm—a sampling strategy— and an L-shaped algorithm to solve the stochastic programming model for a large number of scenarios.
In reality, the cyber network security against attackers can be viewed as a two-player Stackelberg game, where the network owner tries to minimize the loss, whereas the attacker seeks to maximize the loss to the defender. Therefore, to obtain a robust countermeasure deployment decision, the network owner needs to account for the attacker’s objective, making the problem a bi-level optimization. As the network owner has uncertainty about the attackers’ capabilities, some attackers having large capabilities (e.g., resources) can cause severe cyber-attacks, resulting in substantially large financial and reputation loss. Risk-neutral optimization models ignore the risk of these severe attacks. But, these severe cyber-attacks can cripple an organization. Because it is impossible to recover from such consequences in a short time; if recovery is possible at all. Therefore, it is crucial to account for the minimization of the risks of substantially large losses in the presence of uncertainty through a risk-averse modeling approach.
To provide decision-support to the risk-averse network managers, we developed a new risk-averse bi-level stochastic integer programming model that explicitly models (accounts for) multiple attackers with uncertainty in attackers’ capabilities and the risks of substantially large losses from severe attacks while computing the optimal decision [2]. In the stochastic programming model, we incorporated a risk-measure, conditional-value-at-risk (Figure 6), to model the risk of substantially large losses. We developed a customized constraint-and-column generation algorithm as well as several novel speed-up techniques in solving the risk-averse bi-level model. The proposed method provides optimal countermeasure deployment decisions that minimize the expected loss as well as the risk of substantially large losses from severe attacks. Results show that the stochastic risk-averse model provides substantially better network interdiction (countermeasure deployment) decisions (yielding up to 44.78% less expected loss from cyber-attacks) than the existing models that ignore uncertainty in attackers’ capabilities and the risk of severe attacks.
These countermeasure deployment decisions obtained from the optimization methods can then be mapped back in the cyber-network topology to determine on which network devices or links of an organization a security countermeasure should be deployed to protect critical assets. These methods can be easily adapted to secure different smart systems, such as smart manufacturing and smart power grids, from cyber-physical attacks.
Figure 7: Effect of attackers' budget distribution on risk measure
Detection and Simulation of Malicious Cyber-Networks (Botnet) in Large-Scale Networks
Besides the mathematical optimization methods, I worked on developing/implementing supervised machine learning algorithms to detect malicious cyber nodes in large (millions of nodes graph), real-life, and highly class-imbalanced networks utilizing high-performance computing [3]. Malicious entities in the internet network, commonly known as bots, are a serious threat to internet network security, as these bots can perform a wide range of malicious activities such as sending spam emails, conducting distributed denial of service (DDoS), and click-fraud scams. It is challenging to detect bots in large, highly class-imbalanced network data, where the bot proportion is very small compared to the normal traffic. Most of the existing machine learning algorithms fail to detect bots (minority class) in highly class-imbalanced data, as the algorithms are biased towards the majority class. Moreover, most of the existing methods can only detect specific types of bots. Various types of bots can exist in a network, necessitating a generalized bot detection method capable of detecting various types of bots.
Therefore, we developed a novel bot detection method for large (millions of nodes), highly class-imbalanced (bot proportion is very low: 0.00583%), real data. We preprocessed the data to remove noisy observations. By analyzing the data, we developed and computed three novel features that can distinguish bots from normal traffic. Using these features, we developed the following supervised machine learning algorithms: Quadratic discriminant analysis, Gaussian Naïve Bayes, Support Vector Machine, K-Nearest Neighbor, and Random Forest. We trained these algorithms and tuned their parameters to reduce their bias towards the majority class. Due to large-scale (millions of node graphs) data, we used parallel computing to implement our methods in a high-performance computing environment. Results show that the proposed method substantially outperforms the existing methods in detecting bots from large, real-world, highly class-imbalanced data.
In another work, we developed a scalable distributed graph simulation algorithm using the concept of role mining, transition probability matrix, and parallel computing to generate very large-scale graphs (in the scale of billions of nodes) containing malicious nodes [4]. This distributed algorithm is capable of generating very large-scale graphs with accurate simulation of the behavior of malicious nodes from a given real-life input network. Figures 8, 9, and 10 demonstrate the original, simulated, and scaled-up simulated graphs, respectively. We also developed a hypergraph generation algorithm for a parallel computing environment to generate large hypergraphs [5].
We further expand our research on cybersecurity to cyber-physical systems (CPS), which integrate information technology (IT), the Internet of Things (IoT), and operational technology (OT). While this integration has significantly improved the productivity and efficiency of smart manufacturing systems, it has also increased their exposure to cyber-physical attacks. Traditional risk-assessment methods are often subjective, while existing attack-graph models typically assume perfectly effective countermeasures and similar attacker capabilities. To address these limitations, we study the optimal allocation of defense countermeasures in a resource-constrained smart manufacturing system. Our case study is based on the Amatrol 870-Series Mechatronics (AM-870) system at the Cybersecurity for Manufacturing (C4M) Hub, as shown in Figure 1. We thank Mr. Alvaro Espinosa and Mr. Candan Martin of the Cybersecurity of Manufacturing Innovation Institute (CyManII) for their guidance and valuable discussions on the AM-870 case study.
Figure 1. AM-870 case study at the Cybersecurity for Manufacturing (C4M) Hub
Our work is the first study to jointly consider imperfect and multi-level countermeasures, uncertain post-defense attack success likelihoods, and heterogeneous attackers in smart manufacturing. The defender selects countermeasures with different costs and effectiveness levels to minimize the expected loss under a limited defense budget. Figure 2 shows the AM-870 attack graph, which captures the potential paths available to attackers seeking to compromise critical system components.
Figure 2: AM-870 attack graph.
We formulate the problem as a two-stage stochastic bi-level mixed-integer nonlinear programming model based on a single-leader, multi-follower Stackelberg game. The model captures uncertainty in attack success likelihoods and attacker budgets. We reformulate it as a stochastic bi-level mixed-integer linear programming model using probability-chain reformulation and McCormick linearization. We then develop an exact column-and-constraint generation algorithm and a parallel hybrid CCG-PHA-GAA algorithm for solving realistic-sized instances.
Computational results show that the parallel hybrid algorithm is 73 times faster than the column-and-constraint generation method while maintaining a 0.01% optimality gap. Experiments on the AM-870 smart manufacturing system demonstrate that modeling imperfect, multi-level countermeasures produces better defense decisions than assuming binary and perfectly effective defenses. The proposed stochastic model reduces expected loss by up to 41.82% compared with deterministic models. Increasing the number of countermeasure levels and attack-success likelihood states further reduces expected loss by up to 6.28% and 8.96%, respectively.
Bhuiyan, T. H., Nandi, A. K., Medal, H., and Halappanavar, M. (2016). Minimizing expected maximum risk from cyber-attacks with probabilistic attack success. In IEEE Symposium on Technologies for Homeland Security, pp. 1–6. [Available here]
Bhuiyan, T. H., Medal, H., Nandi, A. K., and Halappanavar, M. (2021). Risk-averse bi-level stochastic network interdiction model for cyber-security risk management. International Journal of Critical Infrastructure Protection, 32: 100408. [Available here]
Harun, S., Bhuiyan, T. H., Zhang, S., Medal, H., and Bian, L. (2018). Bot classification for real-life highly class-imbalanced dataset. In 15th IEEE International Conference on Dependable, Autonomic and Secure Computing, pp. 565–572. [Available here]
Bhuiyan, T. H., Harun, S., Medal, H., Zhang, S., (2022). Scalable simulation algorithm to generate large-scale botnet dataset using role-mining and Markov chain. Computers & Security (Under review). [Available here]
Jenkins, L., Bhuiyan, T. H., Harun, S., Lightsey, C., Mentgen, D., Aksoy, S., Stavcnger, T., Zalewski, M., Medal, H., and Joslyn, C (2018). Chapel hypergraph library (CHGL). In IEEE High-Performance Extreme Computing Conference, pp. 1–6. [Available here]
Nguyen, D.N., Bhuiyan, T.H. and Ciocarlie, G. (2026). A Parallel Stochastic Bi-Level Method for Proactive Cyber-Risk Mitigation in Smart Manufacturing Under Endogenous and Exogenous Uncertainties. (Under review).
The rapid growth in e-commerce and demand for timely delivery of perishable or time-sensitive products have underscored the need for faster last mile delivery solutions. Aerial drones offer a distinct potential to improve the last mile delivery of time-sensitive products (e.g., pharmaceuticals, blood samples, and prepared food). However, despite having high speed and energy-efficiency, drones suffer from limited battery capacity and package weight carrying capacity that limit their operability in serving customers far from the distribution center (i.e., depot). Therefore, it is crucial to develop efficient solutions for drone deployment in last mile delivery of time-sensitive products, accounting for the power consumption and battery capacity of drones, requiring battery replacement to serve geographically distributed customers.
Drone Fleet Deployment and Scheduling for Time-Sensitive Deliveries
In one of the works, we studied a deployment optimization problem for a homogeneous fleet of drones performing direct delivery to deliver packages with distinct release times and package pickup times from a depot to geographically distributed customers [1] as shown in Figure 1. The technological developments in drone industry, offer a wide variety of drone types with distinct characteristics—speed, power consumption, battery capacity, payload capacity, and cost. Therefore, we extended this research in our second work for efficient deployment of a mixed fleet of drones with different drone types (shown in Figure 2) for blood sample delivery [2]. We proposed efficient mixed-integer linear programming models, for routing and scheduling of drones with endogenous battery replacement decisions in an as-needed fashion to ensure safe and efficient operation of drones in serving geographically distributed customers. We also developed three solution approaches, including (1) solving the MILP model, enhanced with novel valid inequalities, with CPLEX; (2) a greedy heuristic algorithm; and (3) a heuristic Genetic algorithm.
In one of the works, we studied a deployment optimization problem for a homogeneous fleet of drones performing direct delivery to deliver packages with distinct release times and package pickup times from a depot to geographically distributed customers [1] as shown in Figure 1. The technological developments in drone industry, offer a wide variety of drone types with distinct characteristics—speed, power consumption, battery capacity, payload capacity, and cost. Therefore, we extended this research in our second work for efficient deployment of a mixed fleet of drones with different drone types (shown in Figure 2) for blood sample delivery [2]. We proposed efficient mixed-integer linear programming models, for routing and scheduling of drones with endogenous battery replacement decisions in an as-needed fashion to ensure safe and efficient operation of drones in serving geographically distributed customers.
We also developed three solution approaches, including (1) solving the MILP model, enhanced with novel valid inequalities, with CPLEX; (2) a greedy heuristic algorithm; and (3) a heuristic Genetic algorithm.
We used actual drone flight test data conducted at the Idaho National Laboratory test site and performed a machine learning method to estimate the energy consumption of different types of drones during different flight segments (shown in Figure 3). Numerical results based on a real-life case study of prepared food delivery data from the San Francisco Bay Area, United States, and a real-life case study of actual blood sample delivery from Pendleton, United States, show that the Greedy algorithm solves the problem considerably faster (less than one seconds) than the other approaches, where CPLEX solver is unable to solve the problem within an hour. Results also demonstrate that the mixed fleet include drones with larger package weight carrying capacity and higher battery capacity as the package weight distribution includes heavier packages and drones need to detour through waypoints and fly over the road networks to avoid no-fly zones compared to flying in a straight path as shown in Figure 4.
Figure 2: Different types of drones.
Integrated Drone–Ground Vehicle Routing for Time-Sensitive Deliveries
Drones have limited delivery range and are limited to carry single package per trip, whereas ground vehicles (GVs) despite being less energy-efficient than drones, have long delivery range and are capable of carrying multiple packages in a single trip. Therefore, in our third study [3], we combined the strengths of drones and GVs and proposed an efficient mixed-integer linear programming model for an unsynchronized routing of a mixed fleet of drrones and a homogeneous fleet of GVs to deliver packages with distinct package weights, package release times, and delivery due times to geographically distributed customers as shown in Figure 5.
To provide decision-support tools for logistics system owners, we developed an exact branch-and-cut (B&C) algorithm and a B&C integrated clustering-based heuristic algorithm to solve this problem. Results show that our exact B&C algorithm is up to 22 times faster than the Gurobi solver and our B&C-based heuristic solves the largest problem sizes efficiently while sacrificing the solution quality by a small amount. Results also demonstrate that using only drones in the fleet increases the cost by 141.4% compared to a mixed fleet of drones and GVs, whereas it reduces the energy consumption by 96%. We also studied the effect of weather conditions (e.g., wind speed, precipitation, and temperature) on the operability of different drone types based on a one-year historical data of weather conditions of the San Francisco Bay Area as shown in Figure 6. Results show that the total cost and total energy consumption increase by 19.6% and 13.2%, respectively, as the weather changes from normal to adverse.
Drone Delivery Network Design for Time-Sensitive Products
Using GVs to serve customers that are not directly reachable by drones limits the efficiency and autonomy of drones in the context of last mile delivery. Intermediate battery swapping machines (ABSMs) offer another viable alternative to extend the delivery range of drones by allowing frequent en-route battery swapping of drones to serve geographically distributed customers from a single distribution center (DC). Therefore, in our fourth study [4], we proposed a delivery network design that allows the mixed use of truck delivery and drone delivery supported by ABSMs to serve customers with uncertain demands as shown in Figure 7. We modeled each ABSM as a semi-open queuing network (shown in Figure 8) and adopted matrix geometric method (MGM) to evaluate the performance of ABSMs under different circumstances.
We developed three solution methods, including a heuristic Genetic algorithm, an exact cutting plane algorithm, and a step-wise linear approximation MILP model to solve the delivery network design problem. The numerical results demonstrate that solving the linearized MILP model with Gurobi results in 11.4 times faster runtime compared to the cutting-plane algorithm. Moreover, our proposed customized GA provides high-quality solutions and is 3.6 and 14.3 times faster than the linearized MILP model and the cutting-plane algorithm, respectively. Results also demonstrate that the cost-efficiency of different drone types varies across different demand rates as shown in Figure 9; using the shortest-range drones in the fleet is 39% more cost-efficient than the longest-range drones for the lowest demand rate whereas using the longest-range drones in the fleet is 61% more cost-efficient than the shortest-range drones for the highest demand rate.
Bhuiyan, T. H., Walker, V., Roni, M., and Ahmed, I. (2024), Aerial drone fleet deployment optimization with endogenous battery replacements for direct delivery of time-sensitive products. Expert Systems with Applications, 252 (2024): 124172. [Available here]
Bhuiyan, T.H., Dolatabadi, S.H.H., Uddin, M.J. and Walker, V., 2026. Optimization of a Mixed Fleet of Aerial Drones for Medical Supplies: A Case Study of Blood Delivery Logistics. Transportation Research Record, 2680(8), pp.455-481. [Available here]
Dolatabadi*, S.H.H., Bhuiyan, T.H., Kaleem, W., Subramanyam, A. and Cokyasar, T., 2026. A branch-and-cut algorithm for routing a heterogeneous drone-ground vehicle fleet to deliver time-sensitive products. Transportation Research Part C: Emerging Technologies, 183, p.105454. [Available here]
Dolatabadi, S. H. H., Chowdhury, R., Bhuiyan, T. H., and Dong, W., 2027. Data-driven optimization for drone-truck delivery network design under stochastic demand and semi-open queuing network. Expert Systems with Applications, 334, 134577. [Available here]
Dolatabadi, S.H.H., Bhuiyan, T.H., Zeng, B., Kennedy, J.T. and Neill, B.O., 2026. Designing Hierarchical Hub-and-Spoke Drone-Based Networks for Delivering Time-Sensitive Healthcare Items. arXiv preprint arXiv:2607.12198. [Available here]
The increasing frequency and severity of natural disasters have led to widespread power outages, affecting public health and safety as shown in Figure 1. Therefore, rapid restoration of critical facilities around disasters is vital. However, traditional methods, such as network reconfiguration and dispatching repair crews, are slow and labor-intensive. Moreover, the limited battery capacity of electric vehicles (EVs) and widely varying availability and readiness of personal EVs and electric buses (EBs) necessitate a faster and more reliable solution for restoring power to critical loads, which are disconnected from the grid during a disaster. In our first study [1], we investigate the coordinated routing and energy scheduling of electric school buses (ESBs) to provide temporary power to critical facilities while conventional restoration efforts are underway.
Electric School Bus Routing and Scheduling for Disaster Power Restoration
This is the first study to leverage efficient routing and energy scheduling of a heterogeneous fleet of electric school buses (ESBs) with their robust design, large battery capacity, and structured availability as government-owned assets for faster power restoration of critical isolated loads around disasters as shown in Figure 2. We addressed the following practical aspects: combined transportation and energy scheduling of ESBs, continuous tracking of the ESB’s state of charge (SOC), multiple back-and-forth trips of ESBs between isolated loads (referred to as "shelters") and charging stations (CSs), and route interdependency and spatial-wise coupling among multiple ESB routes. The goal of the decision-maker in this problem is to efficiently route and schedule the ESB fleet that minimizes the total cost including the total cost of the ESBs utilized, and the total transportation energy consumption cost while satisfying the total demand of all the shelters.
We formulated the proposed restoration problem as a novel efficient MILP model with routing and energy scheduling decisions of ESBs. To provide fast and accurate decision-support tools for emergency management agencies, we developed an efficient exact branch-and-price (B&P) algorithm and a customized heuristic B&P algorithm integrating dynamic programming and labeling algorithms.
The numerical results based on a realistic case study of San Antonio disaster shelters and critical facilities demonstrate that our proposed exact and heuristic B&P algorithms are 121 times and 335 times faster, respectively, than the Gurobi solver. Our proposed B&P-based solution algorithms better exploit the network sparsity feature (i.e., limited shelter-ESB type compatibility due to technology, infrastructure, or spatial limitations) than the Gurobi-based approach, enabling faster computations and handling larger problem instances across varying network sparsity levels.
Furthermore, the effective usable capacity metric indicates that the largest-capacity ESBs are 7 times more effective than the least-capacity ESBs. Results also demonstrate that weather severity significantly affects the required fleet size; the required fleet size increases by 225% as the weather changes from normal to adverse. Moreover, our findings show that as the inconvenience fee of demand shifting increases from no inconvenience fee ($0/kWh) to an extremely high inconvenience fee ($10,000/kWh), the percentage of shifted demand reduces to zero, whereas the total restoration cost and the required fleet size increase by up to 200% as shown in Figure 3.
In our next study [2], we extend our earlier ESB routing and energy scheduling research to address uncertainty in shelters’ energy demand and the failure of a single ESB during disaster response. As illustrated in Figure 4, ESBs travel from a depot to serve isolated shelters, make repeated trips to a charging station, and return to the depot before the planning horizon ends. While the earlier study determines a heterogeneous fleet and its operations to satisfy known energy demand at minimum restoration cost, this study considers a fixed, available fleet of homogeneous ESBs and recognizes that limited resources and disruptions may prevent full demand satisfaction. We therefore seek to minimize the inconvenience cost of unmet demand under worst-case conditions, while retaining continuous battery state-of-charge tracking and coordination among ESBs serving shared shelters.
We formulated the problem as a two-stage robust optimization model, where ESB routes and service schedules are determined before uncertainties are observed, while energy delivery decisions adapt to realized shelter demand and single-ESB failures. To solve the model efficiently, we developed an exact decomposition algorithm combining column-and-constraint generation (C&CG) and Lagrangian relaxation (LR). We further accelerated the approach by introducing additional decomposition layers that divide the Lagrangian subproblems into smaller subproblems with analytical solutions, providing a faster exact decision-support tool for emergency management agencies.
Numerical results based on a real-life case study of disaster shelters in San Antonio, Texas, demonstrate that the proposed Lagrangian-based C&CG algorithm with analytical subproblem solutions is 672 times faster than basic C&CG. The results also reveal a non-monotonic relationship between shelter service time and worst-case unmet demand: increasing service time from one to three time slots reduces the percentage of worst-case unmet demand by 63%, whereas increasing it from three to five time slots raises this percentage by 457%. Furthermore, increasing the available fleet from two to six ESBs reduces the percentage of worst-case unmet demand by 65%, highlighting the importance of fleet availability and service-time planning for emergency power delivery under uncertainty.
Proactive disaster preparedness and resource allocation for resilient power distribution systems
In our third study [3], we broaden our focus to proactive preparedness and reactive restoration of power distribution systems under extreme weather events. From the perspective of a distribution system operator, we develop a two-stage distributionally robust optimization framework that coordinates the allocation of limited resilience resources before disasters with network restoration during disasters. The first stage captures multi-period preparedness decisions, including line hardening, while the second stage coordinates charging station schedules to support critical loads and operate the damaged network. To address uncertain electricity demand and infrastructure contingencies, we introduce a joint decision-dependent ambiguity set that represents a range of plausible probability distributions and captures how preparedness decisions affect network vulnerability. This coupled framework links pre-disaster investments with subsequent restoration performance, enabling the operator to plan against adverse probability distributions without relying on precisely known uncertainty distributions.
Dolatabadi, S.H.H., Dong, Y., Bhuiyan, T. H., Zeng, B., ONeill, B. and Severson, A. (2025). Leveraging Electric School Buses for Disaster Recovery: Optimizing Routing and Energy Scheduling via Branch-and-Price. arXiv preprint arXiv:2510.14131. [Available here]
Dolatabadi, S.H.H., Bhuiyan, T. H., Zeng, B. (2026). A Nested Lagrangian-Based Column-and-Constraint Generation Algorithm for Disaster Restoration Leveraging Electric School Buses. Working paper
Nguyen, D.N., Bhuiyan, T. H. and Zeng, B. (2026). Proactive disaster preparedness and resource allocation for resilient power distribution systems. Working paper
The transition toward decentralized and integrated energy systems creates new challenges for coordinating energy transactions among diverse market participants. Effective market mechanisms are therefore needed to coordinate distributed energy resources while balancing the economic and operational interests of different stakeholders. This motivates the development of advanced market designs and decision-support frameworks for the efficient and sustainable operation of increasingly complex local energy systems.
Local Electricity Market Design and Management
The growing adoption of electric vehicles (EVs) and renewable energy sources creates opportunities for cleaner and more flexible electricity systems. However, managing local electricity markets (LEMs) is challenging because participants have different resources and financial interests, while electricity generation and consumption are uncertain. In our first study [1], we developed a stochastic game-theoretic optimization framework for managing a competitive LEM that integrates EVs, renewable generation, and energy storage, as shown in Figure 1. A local transaction center (LTC) acts as the market operator, setting distinct hourly electricity prices for a load aggregator representing residential complexes and a charging station managing solar generation, battery storage, and EV charging. Each entity seeks to maximize its own profit while satisfying its operational requirements.
We formulated the problem as a bi-level, two-stage stochastic optimization model incorporating heterogeneous participant behaviors, demand response, vehicle-to-grid participation, and uncertainty in generation and consumption. Day-ahead prices and energy transaction commitments are determined before uncertainty is realized, while operational decisions adapt to realized conditions. We developed a centralized solution approach based on Karush–Kuhn–Tucker (KKT) conditions and a decentralized approach using a hybrid genetic algorithm. The centralized approach provides optimal solutions, while the decentralized approach reduces computation time and preserves participants’ private information by allowing them to solve their own operational problems.
Numerical results from a realistic case study demonstrate how dynamic pricing coordinates local energy trading. As shown in Figure 2, participants switch between buying and selling electricity throughout the day as their energy needs and available resources change. Figure 3 shows that the proposed dynamic pricing mechanism enables profitable operation for all three market entities in the tested setting, whereas unified hourly and flat pricing result in losses for the LTC. Figure 4 compares participants’ profits across the solution approaches, with the exact KKT approach achieving the highest LTC profit. For the 10-scenario case, the hybrid genetic algorithm reduces computation time from approximately 3.8 hours to 44 minutes, with a 21% optimality gap, demonstrating the trade-off between computational efficiency and solution quality.
Figure 4: comparison of market participants’ profits across the solution approaches.
Further experiments highlight the benefits of energy storage and uncertainty-aware planning. Increasing the LTC’s storage capacity from 200 MWh to 400 MWh, while changing its initial condition from fully discharged to fully charged, increases its operating profit by 255% and reduces dependence on the main grid without changing the other participants’ profits. At the highest uncertainty level tested, stochastic planning produces 433% higher expected LTC profit than decisions based on average forecasts, when evaluated under uncertain conditions with penalties for energy imbalances. These findings demonstrate the importance of dynamic pricing, flexible storage resources, and explicit uncertainty modeling in designing financially viable local electricity markets.
Parallel Decentralized Algorithm for Coupled Electricity and Thermal Market Management
We extend the electricity-only market framework to a hierarchical local coupled power and heat market (H-LCPHM), addressing the coordinated operation of electricity and thermal resources under increasing renewable energy integration [2]. As illustrated in Figure 5, the expanded market incorporates combined heat and power units and heat pumps, coupling electricity and heat to improve system efficiency, flexibility, and sustainability. The framework captures heterogeneous market participants, including a local transaction center (LTC), load aggregators (LAs), industrial prosumers (IPs), and their lower-level prosumers, each seeking to maximize its own profit within a competitive game-theoretic structure. We combine dynamic pricing with coordinated electricity and heat management and develop a Karush–Kuhn–Tucker (KKT)-based decomposition approach to solve the resulting market optimization problem.
Figure 5: Hierarchical local coupled power and heat market structure.
Dolatabadi, S.H.H., Bhuiyan, T. H., Chen, Y. and Morales, J.L., 2024. A stochastic game-theoretic optimization approach for managing local electricity markets with electric vehicles and renewable sources. Applied Energy, 368, p.123518. [Available here]
Nguyen, D.N., Bhuiyan, T. H., Zeng, B. and Chen, Y., 2026. A Parallel Decentralized Algorithm for Managing Hierarchical Coupled Local Electricity and Thermal Markets with Electric Vehicles and Renewable Sources. Working paper
Many scientific and engineering problems require decision-making to learn unknown system behavior or identify optimal solutions. However, in extreme data-scarce environments, conventional decision-making and search methods are inefficient because they always require a large number of expensive sample or function evaluations, resulting in substantial experimental or computational costs. Bayesian active learning (BAL) provides an effective approach for sequential decision-making by determining which samples to evaluate based on the estimated mean and uncertainty of a probabilistic surrogate model. This makes BAL particularly valuable in data-scarce, black-box settings where the underlying functional relationships are unknown, each sample evaluation is expensive, and only a limited number of evaluations can be performed. BAL generally serves two primary goals: optimization (also known as Bayesian optimization), which focuses on efficiently identifying optimal solutions, and surrogate model learning, which focuses on improving the predictive accuracy of the surrogate model over the decision space.
Batch Multi-Objective Bayesian Optimization for Data-Efficient Refractory High-Entropy Alloy Design
In one of our projects, we studied a multi-objective refractory high-entropy alloy (RHEA) design problem for nuclear reactor applications. RHEA design involves a large compositional design space, expensive experimental evaluations, multiple conflicting objectives, and black-box relationships between the design variables and material properties (i.e., closed-form functional relationships are unavailable). As a result, exhaustive expensive experimentation is impractical, and conventional optimization methods that rely on explicit functional relationships cannot be applied. Therefore, in one of our projects, we proposed a batch multi-objective Bayesian optimization (BMOBO) method for the data-efficient discovery of RHEAs under multiple conflicting material properties (i.e., design objectives) [1].
Our proposed BMOBO frameworks can efficiently select multiple promising alloy compositions and reactor operating temperatures while balancing four conflicting design objectives: yield strength, elongation, Vickers hardness, and oxidation resistance. To account for different representations of reactor operating temperatures, we developed two BMOBO frameworks: a continuous-design-variable framework (BMOBO-CV), which treats temperature as continuous, and a mixed-design-variable framework (BMOBO-MV), which treats temperature as ordinal to better reflect practical reactor operating conditions. The proposed frameworks balance exploration and exploitation to efficiently identify improved Pareto-optimal alloy designs under a limited experimental budget.
We first evaluate the frameworks using datasets from the literature. After five optimization iterations, BMOBO-MV outperformed BMOBO-CV, improving proportional hypervolume (PHV) by 18.27% and reducing inverted generational distance (IGD) by 53.85%. Compared with the conventional Non-dominated Sorting Genetic Algorithm, the BMOBO frameworks improved PHV by 32.26% and reduced IGD by 60.00%, demonstrating the advantage of BMOBO-MV in efficiently exploring the RHEA design space under limited data and evaluation budgets. Figure 1 illustrates the improvement achieved by BMOBO-MV from the first BAL iteration to the fifth BAL iteration (i.e., the final BAL iteration), showing how the identified Pareto front progressively improves over the sequential decision-making process.
Figure 1: Pareto front across five BMOBO-MV iterations.
We further conducted experimental validation of our proposed BMOBO-MV framework. At each BAL iteration, the framework recommended an alloy composition and reactor operating temperature, after which the corresponding RHEA sample was fabricated and experimentally characterized. Microstructural and compositional characteristics were examined using scanning electron microscopy (SEM) and energy-dispersive X-ray spectroscopy (EDS), while hardness and oxidation resistance were experimentally measured. Figure 2 shows the microindenter used for sample hardness testing, and Figure 3 shows the EDS results of Sample 1. Experimental validation showed that BMOBO-MV provides Pareto-optimal alloys within only two sequential evaluations that achieve hypervolume and individual alloy properties by up to 211.52% and 75%, respectively.
Cost-Aware Bayesian Active Learning with Derivative Information for Surrogate Modeling
In another project, we studied the problem of accurate surrogate model construction under limited data-acquisition budgets in data-scarce environments (e.g., thermal analysis of a heated fin, additive manufacturing problems) using the proposed cost-aware BAL framework [2]. Unlike optimization-oriented BAL, the goal of this project is not to identify an optimal solution, but to learn a better surrogate model that approximates the unknown function over the entire decision space under a limited data-acquisition budget. In addition to conventional function-value observations, derivative observations can provide complementary information about the local behavior of the unknown function, but these different observation types may incur different evaluation costs.
To leverage the additional information provided by derivative observations, we developed a cost-aware generalized variance reduction (GVR) acquisition strategy for BAL that quantifies the expected reduction in global surrogate-model uncertainty relative to the corresponding evaluation cost of different observation types. By balancing the information gained from each observation against its evaluation cost, the framework aims to achieve a desired level of surrogate-model accuracy with a lower total evaluation cost. Figure 4 compares the effects of different integration resolutions used in the GVR calculation on predictive accuracy, computation time, and total evaluation cost.
Figure 4: Effect of Integration Resolution on Cost-Aware GVR
Graham, K., Li, K., Bhuiyan, T. H., Warren, P., Khair, M., Sooby, E., & Millwater, H. (2026). A batch multi-objective Bayesian optimization approach for nuclear high-entropy alloy design. Working paper.
Li, K., Bhuiyan, T. H., Morales Bernal, D. F., Roberts, S., Restrepo, D., & Millwater, H. (2027). Bayesian Active Learning for Additive Manufacturing Modeling under Extreme Data Scarcity. Working paper.
In my Ph.D. dissertation, I developed a data-driven SPEU modeling framework where both exogenous and endogenous uncertainties coexist to model an integrated facility protection and transportation network design problem under random disruptions with the following realistic assumptions: (1) effect of protection resources is imperfect, (2) protection is multi-level, and (3) facilities have multiple post-disruption capacities. The proposed two-stage stochastic program integrates supervised machine learning models to estimate the conditional probability distributions of the facilities’ post-disruption capacities by analyzing historical data. To linearize the non-convex, nonlinear stochastic program, we developed a new mathematical reformulation technique. Utilizing the structure of the model, we developed a customized L-shaped decomposition algorithm to solve the problem as well as introduced several new valid inequalities to improve the computational efficiency of the algorithm.
Figure 2: A distribution network subject to disruptions based on the southeastern United States
The numerical results based on a distribution network in the southeastern United States (as shown in Figure 2) demonstrate that the optimal solution from the proposed methods provides an average reduction of 18.71% in expected transportation cost than the existing methods assuming perfect and binary protection and disruption. Also, the proposed methodology provides substantially (on average 600%) lower expected transportation costs than the methods ignoring the endogenous uncertainty in the facilities’ post-disruption capacities. This research provides several key managerial insights about the effect of budget variation on the expected post-disruption transportation cost and the effect of the estimation error of the supervised machine learning algorithms on the optimal solution of the optimization model. These insights are beneficial for the network managers in deciding on the investment amount in reliable network design and critical infrastructure protection as well as in choosing the machine learning algorithms that have the least impact on the optimization model.
Bhuiyan, T. H., Medal, H. R., and Harun, S. (2020). A stochastic programming model with endogenous and exogenous uncertainty for reliable network design under random disruption. European Journal of Operational Research, 285(2):670–694. [Available here]
Centralized government (land and fire management) agencies often face challenges in implementing wildfire prevention measures (hazardous fuel reduction treatments, such as prescribed burning, mechanical or chemical vegetation control) on nonindustrial private forest lands. These agencies can only incentivize the private landowners to encourage them to implement wildfire prevention measures on their lands through a cost-share program by sharing some portion of the implementation cost. The key challenges for the agencies in designing a cost-share program with their limited budget are: (1) the agency does not know whether a landowner will accept or reject the cost-share program for a given financial assistance amount, (2) different landowners’ land have different biophysical characteristics (e.g., vegetation density, moisture content, elevation, slope, etc.) and pose different amount of wildfire risk to the surrounding landscape.
In my Ph.D. dissertation, I proposed a novel risk-based incentive structure design program to incentivize private landowners in implementing fuel reduction treatments on their lands, where (1) the landscape is divided into a number of raster cells (shown in Figure 3) based on the biophysical characteristics, and (2) different landowners are offered different incentive amounts by the budget-constrained government agency depending on the wildfire risks each landowner's land pose to the surrounding landscape. Another challenge the government agency faces is that the agency does not know whether a landowner will accept or reject the cost-share program for a given incentive amount. However, based on our landowner survey data, we found that the landowners are more likely to accept the cost-share program as the incentive amount increases, which makes the landowners' behavioral uncertainty endogenous to the incentive offer decision.
We modeled the risk-based incentive structure design problem in reducing the risk of catastrophic future wildfires as a data-driven SPEU model. The proposed SPEU model accounts for the endogenous uncertainty in the landowners’ accept/reject decisions of the incentive. The proposed SPEU model utilizes a supervised machine learning model to estimate the landowners' likelihood of accepting the cost-share program by analyzing landowner survey data. Due to a large number of scenarios, the second-stage value of the two-stage SPEU model is computed using a simulation program that models the wildfire spread in the landscape accounting for the information on fuel reduction treatment, weather, and landscape characteristics. To model the wildfire spread as a network flow problem, the landscape is represented as a directed network (Figure 4), where each node represents a land parcel of the landscape and the arc weights represent the fire travel time between the land parcels. This work falls in the class of stochastic network interdiction problems, as the spatial optimization of fuel treatment resources prevents future catastrophic wildfires.
Results based on real data from Santa Fe National Forest, New Mexico show that the proposed method can provide decision support to land management agencies for a realistic-sized landscape. This novel risk-based approach provides up to 37.3% more reduction in expected wildfire damage than the currently used methods.
Bhuiyan, T. H., Moseley, M. C., Medal, H. R., Rashidi, E., and Grala, R. K. (2020). "A stochastic programming model with endogenous uncertainty for incentivizing fuel reduction treatment under uncertain landowner behavior". European Journal of Operational Research, 277(2):699–718. [Available here]
I am studying the vulnerability of multi-hop wireless networks to jamming attacks. Gaining insights into network vulnerability is the first step in protecting wireless networks from adversaries. Unlike wired networks (e.g., transportation networks as shown in Figure 13), wireless networks (Figure 14) have a special property called radio interference that significantly affects the network throughput and makes the modeling of wireless networks to be different than wired networks. Due to this radio interference property, communications interfere with each other if they are within their range of interference and thus all communications cannot be active simultaneously, which substantially affects the network throughput. Figure 14 shows a connectivity graph of a small wireless communication network, where the same colored arcs can communicate (be active) simultaneously, which forms a schedulable set, a.k.a., independent set in graph theory. Based on this interference among communications, we can construct a conflict graph (or interference graph) to represent the radio interference. Figure 15 demonstrates a partial conflict graph corresponding to the connectivity graph shown in Figure 14. An independent set of nodes (equivalently a set of arcs in the connectivity graph) in the conflict graph can transmit data simultaneously, and scheduling maximal independent sets ensures maximum flow in the network.
There exist several models of this radio interference, such as capture, protocol, and interference range. However, However, a physical model of radio interference in wireless networks is more realistic; but, the jammer placement problem under the physical model of interference has not been analyzed. We have developed attacker-defender bi-level mixed-integer programming models under the physical (i.e., additive) and other non-additive models (e.g., capture, protocol, interference range) of radio interference to find optimal locations of jamming devices that minimize the maximum network throughput. We have developed an improved branch-and-cut and a hybrid Tabu search algorithm to solve the bi-level models. This research provides key insights and helps in finding the most beneficial strategies for designing a network to be jamming-resistant and can lead to new tools for designing wireless networks to mitigate the effect of jamming attacks.
Bhuiyan, T. H., Medal, H. (2022). Computing the vulnerability of multi-hop wireless networks: A search for the "best" interference model. Working paper.
I applied the SPEU solution approaches developed in my Ph.D. research to mitigate supplier risks—poor quality and delivery lateness—for manufacturing companies. We studied the problem of selecting supplier development programs (SDPs) to proactively mitigate supplier risks, where the effects of SDPs on supplier performance improvement are uncertain and depend probabilistically on the SDP selection decision. We proposed a new SPEU model that selects optimal SDPs accounting for multiple supplier risks [1]. We implemented an accelerated L-shaped decomposition algorithm and proposed a sample-based greedy approximation algorithm to solve the model. While the L-shaped decomposition algorithm provides the exact solution for moderately-sized problems, the sample-based greedy approximation algorithm provides very good-quality approximate solutions for large-scale problem instances. Our case study based on four manufacturing companies shows that the proposed model and solution approaches are beneficial for low-volume, high-value manufacturing industries (e.g., nuclear power plant construction, shipbuilding).
Zhou, R., Bhuiyan, T. H., Medal, H., Sherwin, M. D., and Yang, D. (2022). A stochastic programming model with endogenous uncertainty for selecting supplier development programs to proactively mitigate supplier risk. Omega, 107: 102542. [Available here]