Stochastic Optimization of Solvent Recovery Networks with Equipment Degradation Dynamics: A Two-Stage Programming Framework for Sustainable Industrial Operations

Abstract

Industrial solvent recovery systems operate at the intersection of economic pressure, regulatory constraint, and physical deterioration—a combination that makes purely deterministic optimization frameworks inadequate for real-world deployment. This article presents a comprehensive methodology for optimizing solvent recovery networks under simultaneous uncertainty in feed composition, market prices, and equipment performance, with explicit attention to the temporal dynamics of equipment degradation. The proposed framework integrates a two-stage stochastic programming model with a Markov-modulated degradation process, enabling both strategic maintenance scheduling and real-time operational recourse decisions. Scenario generation is accomplished through Latin hypercube sampling coupled with moment-matching, and the resulting large-scale stochastic programs are solved via an enhanced L-shaped decomposition method with strengthened Benders cuts. Extensive computational results and two industrial case studies—one from the automotive coating sector and one from pharmaceutical active ingredient manufacturing—demonstrate that the proposed approach yields cost savings of 12–19% relative to nominal deterministic models while maintaining solvent recovery rates above regulatory thresholds even under adversarial realizations of uncertainty. The methodology advances the state of the art by treating equipment health as a first-class decision variable rather than a fixed parameter, thereby connecting operational and maintenance planning within a single coherent optimization model.

Keywords: solvent recovery, process optimization, stochastic programming, equipment degradation, sustainable manufacturing, maintenance scheduling, chemical process systems

1. Introduction

Volatile organic compounds (VOCs) emitted during industrial coating, cleaning, and synthesis operations represent one of the most consequential environmental challenges in modern manufacturing. Solvent recovery—the recapture and reuse of these compounds through condensation, adsorption, or membrane separation—has grown from a compliance-driven afterthought into a core element of sustainable manufacturing strategy. As regulatory frameworks such as the European Solvents Emissions Directive and the U.S. National Emission Standards for Hazardous Air Pollutants continue to tighten, process owners face mounting pressure to achieve high recovery rates while simultaneously managing the economic viability of recovery operations (Capello et al., 2007; U.S. Environmental Protection Agency, 2020).

The economics of solvent recovery are inherently volatile. Recovered solvent value tracks commodity markets for raw materials such as acetone, ethyl acetate, and isopropyl alcohol; these markets exhibit seasonal cycles, supply shocks, and demand fluctuations that can shift break-even recovery rates by tens of percentage points over the course of a planning horizon (Henderson et al., 2011). At the same time, the physical performance of recovery equipment—adsorption beds, condensers, distillation columns, and membrane modules—degrades continuously with use. Activated carbon beds lose capacity as micropores foul with high-molecular-weight contaminants; heat exchangers accumulate fouling films; membrane modules experience compaction and plasticization. Ignoring this degradation in optimization models leads to optimistic performance predictions that fail to materialize in practice, often resulting in unplanned shutdowns and missed regulatory targets (Scarf, 1997).

Despite decades of progress in process systems engineering (PSE), the joint treatment of operational uncertainty and equipment degradation within solvent recovery network optimization remains underdeveloped. Early work by El-Halwagi (1997) established rigorous mass-integration frameworks for solvent recovery, and subsequent extensions addressed multicomponent separation sequencing and network synthesis (El-Halwagi, 2012). Stochastic programming has been applied to chemical process networks in contexts ranging from supply chain management (Gupta & Maranas, 2003) to multisite production planning (Verderame & Floudas, 2010), but degradation dynamics have rarely been incorporated explicitly. In parallel, the reliability engineering and maintenance optimization communities have developed sophisticated degradation models (Nakagawa, 2005; Zio, 2009), but these typically treat the process operating point as fixed rather than as a decision variable.

The present work bridges these traditions by formulating solvent recovery network optimization as a two-stage stochastic program in which: (1) first-stage decisions encompass network configuration, maintenance scheduling, and capacity allocation; and (2) second-stage recourse decisions cover real-time flow routing, auxiliary regeneration, and blending adjustments conditional on observed uncertainty realizations. Equipment degradation is modeled as a discrete-state Markov chain whose transition probabilities depend on operating conditions, directly coupling operational choices to long-run equipment health trajectories. This formulation enables the optimizer to trade off short-term recovery revenue against long-term equipment health in a principled way.

Two industrial case studies ground the methodology in practice. The first examines a multi-stream solvent recovery system in an automotive OEM coating facility, where the primary solvents—xylene, methyl ethyl ketone (MEK), and n-butanol—arrive in highly variable feed concentrations driven by paint formulation changes. The second studies a pharmaceutical synthesis plant where acetonitrile and dichloromethane recovery must satisfy both economic objectives and stringent purity constraints tied to active pharmaceutical ingredient (API) quality. Together, these cases illustrate how the framework adapts to different uncertainty profiles, regulatory environments, and equipment portfolios.

The remainder of the article is organized as follows. Section 2 provides a concise review of relevant prior work. Section 3 presents the full mathematical formulation, including network topology, degradation dynamics, stochastic programming structure, and the risk measure used to balance expected performance against tail risk. Section 4 describes the solution algorithm and scenario generation procedure. Section 5 presents computational results and case study analyses. Section 6 discusses broader implications and limitations. Section 7 concludes.

2. Background and Related Work

2.1 Solvent Recovery Network Optimization

The systematic design of solvent recovery networks draws heavily from mass integration, a branch of PSE concerned with the optimal routing and reuse of material streams within a process complex (El-Halwagi, 1997). The source–sink matching framework, first popularized in analogy with heat integration (Linnhoff et al., 1982), identifies recovery opportunities by comparing the quality attributes—most importantly, solvent concentration—of available streams against the requirements of potential sinks. El-Halwagi and Manousiouthakis (1989) introduced the concept of mass exchange networks, which formalize the problem as a mathematical program with logical feasibility constraints.

More recently, researchers have extended these ideas to networks involving multiple separation technologies operating in parallel or in series. Seuranen et al. (2005) examined the sequencing of distillation and adsorption units for multicomponent solvent recovery, demonstrating that the choice of separation sequence has first-order effects on energy consumption and recovered product purity. Capello et al. (2007) offered a lifecycle perspective, showing that the environmental credits from solvent recovery can outweigh the energy costs of recovery operations across a wide range of conditions, reinforcing the sustainability rationale for investment in sophisticated recovery infrastructure.

2.2 Stochastic Programming in Chemical Process Systems

Stochastic programming provides a rigorous mathematical language for optimization under uncertainty (Birge & Louveaux, 2011; Shapiro et al., 2014). In the PSE literature, two-stage stochastic programs have been applied to process design under uncertain demand (Acevedo & Pistikopoulos, 1998), production scheduling under processing time uncertainty (Gupta & Maranas, 2003), and supply chain network design under price volatility (Gebreslassie et al., 2012). Ruszczyński and Shapiro (2003) provide the theoretical foundations underpinning the decomposition algorithms commonly used to solve these large-scale problems.

Risk measures have become increasingly important in process industry applications. The conditional value-at-risk (CVaR), introduced by Rockafellar and Uryasev (2000) and applied to process systems by Barbosa-Póvoa (2012), penalizes the tail of the cost distribution rather than only its expectation, making the optimization result more robust to extreme scenarios—a property of clear value when regulatory violations carry disproportionate penalties.

2.3 Equipment Degradation Modeling

Degradation models in engineering range from physics-based mechanistic formulations to data-driven statistical approaches. Nakagawa (2005) surveys classical maintenance models built on exponential and Weibull failure-time distributions. More nuanced continuous-time degradation processes include gamma processes, inverse Gaussian processes, and Wiener processes with drift, each suited to different failure mechanisms (Zio, 2009). For adsorption-based solvent recovery, the degradation of activated carbon capacity has been described empirically by first-order kinetic models linking fouling rate to inlet contaminant concentration and temperature (Ruthven, 1984).

The Markov chain formulation adopted here follows a tradition in infrastructure management (Madanat et al., 1995) and has seen some adoption in chemical plant maintenance planning (Scarf, 1997). Its primary advantage is computational tractability: the transition probability structure lends itself naturally to integration within a mathematical programming framework, avoiding the simulation overhead required by continuous degradation models.

2.4 Joint Operational and Maintenance Optimization

The integration of maintenance scheduling with operational optimization has received growing attention in the PSE community. Pistikopoulos et al. (2000) examined simultaneous production and maintenance scheduling in multipurpose batch plants, showing that ignoring maintenance in scheduling leads to systematic underperformance. Dedopoulos and Shah (1995) addressed the reliability of multiproduct plants under maintenance constraints. More recently, Cabrera-Ríos et al. (2022) applied stochastic dynamic programming to integrated maintenance and operation of distillation systems under demand uncertainty. The present work extends this line of research to recovery networks with heterogeneous equipment portfolios and explicitly stochastic uncertainty structures.

3. Mathematical Formulation

3.1 Network Topology and Mass Balance

Consider a solvent recovery network consisting of a set of source streams , a set of separation and recovery units , and a set of sink streams . Each source stream is characterized by a volumetric flow rate (L/h) and a composition vector where denotes the concentration (g/L) of solvent species . The overall mass balance for species across unit is given by:

where (L/h) is the flow from source to unit , (L/h) is the flow from unit to sink , is the outlet concentration of species directed to sink , and is the loss (unrecovered emission) of species from unit . Equation (1) enforces conservation of mass for each species at each unit, forming the backbone of the network feasibility constraints.

The recovery efficiency of unit for species defines the fraction of incoming species that is captured and directed to recovery sinks. Critically, this efficiency is not a fixed parameter but a function of the unit's current degradation state , where represents a new or fully regenerated unit and represents end-of-life condition:

where is the peak efficiency of the unit (achieved at ) and is a species-specific degradation sensitivity coefficient. This parameterization reflects the empirical observation that high-molecular-weight contaminants disproportionately foul adsorption sites relevant to lighter solvents (Ruthven, 1984). Equation (2) is deliberately kept linear in to preserve the mixed-integer linear programming (MILP) structure of the overall problem; more complex nonlinear relationships can be linearized via piecewise linear approximation if warranted by data availability.

3.2 Equipment Degradation Dynamics

Equipment degradation is modeled as a discrete-time, finite-state Markov chain with state space and time index . The transition probability from state to state over one time period is denoted and depends on the operating intensity —a continuous variable representing the throughput fraction of unit at time . Specifically:

for , while state is absorbing until a maintenance action is taken. The degradation rate is a unit-specific parameter estimated from historical operational data using maximum likelihood estimation. This formulation captures the intuition that higher throughput accelerates wear—a well-documented phenomenon in adsorption and membrane filtration systems (Zio, 2009).

Maintenance decisions are represented by the binary variable , where indicates that unit undergoes maintenance (regeneration or overhaul) at the start of period . When maintenance is performed, the degradation state resets:

where is the state resulting from the Markov transition in the absence of maintenance, and is the post-maintenance state (typically zero for full overhaul, or some residual value for partial regeneration). To maintain linearity, the expected degradation state is used in the optimization: , where is the probability of being in state at time , updated according to the Chapman–Kolmogorov equations.

3.3 Uncertainty Representation and Scenario Generation

Three principal sources of uncertainty affect solvent recovery network operations: (1) feed stream compositions , which fluctuate with upstream production schedules and raw material variations; (2) recovered solvent prices , which track commodity markets; and (3) utility costs (primarily energy for distillation and refrigeration). These are collected into a random vector defined on a probability space .

Scenario generation follows a two-step procedure. First, marginal distributions for each uncertain parameter are fitted from historical process data using kernel density estimation. Second, a Latin hypercube sampling (LHS) scheme with moment-matching correction (Høyland & Wallace, 2001) generates a finite scenario set with associated probabilities summing to unity. The moment-matching step ensures that the discrete scenario tree reproduces the mean, variance, skewness, and correlations of the underlying continuous distributions to within a specified tolerance, which is important for accurately capturing the right tail of the cost distribution—the region most relevant to regulatory risk.

3.4 Two-Stage Stochastic Programming Formulation

The two-stage stochastic program is organized as follows. First-stage decisions (here-and-now, made before uncertainty is resolved) include: network configuration binary variables indicating whether unit is active in the planning horizon; maintenance schedules ; and nominal capacity allocations . Second-stage decisions (wait-and-see, contingent on scenario ) include flow rates and auxiliary operational adjustments .

The objective is to minimize the sum of first-stage costs and the risk-adjusted expectation of second-stage costs:

where denotes first-stage costs (capital investment amortization, scheduled maintenance costs, and fixed operating costs), is the second-stage value function representing net operating cost in scenario , is a risk-aversion weight, and is the conditional value-at-risk at confidence level (typically in our formulations). Setting recovers the risk-neutral (expected value) objective; increasing shifts emphasis toward protecting against worst-case scenarios.

The CVaR is linearized using the standard auxiliary variable formulation of Rockafellar and Uryasev (2000):

where is the value-at-risk threshold variable, is replaced by the auxiliary variable satisfying , yielding a linear reformulation amenable to MILP solvers.

The second-stage value function for scenario is:

where is the mass of solvent recovered in scenario , is the scenario-dependent operating cost (energy, consumables), and is a penalty term for regulatory violations—activated whenever the recovery rate for any regulated species falls below the mandated threshold . The penalty is modeled as a piecewise linear function of the shortfall, reflecting escalating regulatory fines.

3.5 Key Constraints

The complete model includes the following constraint classes:

Flow feasibility and capacity constraints:

Degradation state evolution (mean-field approximation):

Recovery efficiency coupling:

Maintenance resource constraints: At most units can undergo maintenance simultaneously due to crew availability:

Regulatory compliance constraint:

where is the regulatory shortfall that enters the penalty term in Equation (7).

4. Solution Methodology

4.1 L-Shaped Decomposition with Strengthened Cuts

Direct solution of the extensive form of the stochastic program is computationally prohibitive for large scenario sets and long planning horizons. With scenarios, the extensive form contains on the order of binary variables and continuous variables for typical network sizes, well beyond the practical reach of branch-and-bound MILP solvers without decomposition.

The L-shaped method (Van Slyke & Wets, 1969), generalized to problems with integer first-stage variables (Laporte & Louveaux, 1993), decomposes the problem into a master program over first-stage decisions and independent second-stage subproblems. The master program iteratively receives optimality cuts (Benders cuts) from the subproblems:

where superscript denotes iteration number and is the epigraphical variable approximating the recourse cost in the master.

Three enhancements are incorporated to accelerate convergence:

1. Magnanti–Wong strengthening: Subproblem dual solutions are projected to the relative interior of the dual feasible set before cut generation, producing cuts that are active at a richer set of first-stage solutions and typically cutting off larger portions of the feasible space per iteration (Magnanti & Wong, 1981).
2. Multi-cut aggregation: Rather than aggregating all scenario cuts into a single expected-value cut, we maintain separate cuts, which provides a tighter lower bound approximation at the cost of a larger master program. Empirically, this trade-off favors the multi-cut approach for our problem sizes.
3. Warm-starting from deterministic solution: The deterministic nominal problem (using mean values of ) is solved first, and its solution is used to initialize the first-stage variable values in the stochastic master, substantially reducing the number of L-shaped iterations required.

4.2 Rolling Horizon Implementation for Real-Time Operations

While the two-stage stochastic program is well-suited for strategic planning (monthly or quarterly horizons), real-time operational adjustments require a receding horizon implementation. At each operational period , the observed realization is used to update the scenario distribution (Bayesian update of distribution parameters) and a reduced-horizon stochastic program over the remaining periods is solved. This rolling horizon approach (Chand et al., 2002) allows the model to continuously incorporate new information about equipment state—measured via inline sensors tracking adsorption breakthrough curves, heat exchanger pressure drops, or membrane flux decline—while maintaining consistency with the long-horizon maintenance schedule.

The degradation state estimate is updated using a Kalman-filter-inspired correction step that combines the Markov chain prediction with sensor observations (e.g., observed recovery efficiency or pressure drop indicator):

where is the Markov-chain predicted state, maps the observable to a degradation state estimate (derived from the inverse of the efficiency-degradation relationship in Equation (2)), and is a gain parameter that weights model prediction against sensor data, calibrated during a commissioning phase.

4.3 Computational Implementation

The optimization framework is implemented in Python 3.10 using the Pyomo algebraic modeling language (Hart et al., 2011) interfaced with the Gurobi 10.0 MILP solver (Gurobi Optimization, 2023). The L-shaped decomposition is implemented as a custom callback within Gurobi's lazy constraint framework. Scenario generation uses the SciPy LHS implementation for initial sampling and a custom moment-matching routine based on the algorithm of Høyland and Wallace (2001). All computations are performed on a server with 64 Intel Xeon cores and 256 GB RAM, with subproblem solves parallelized across cores using Python's multiprocessing module.

# Illustrative pseudocode for the L-shaped decomposition master iteration
for iteration in range(max_iter):
    # Solve master program
    master_result = solve_master(master_model, solver='gurobi')
    y_k = master_result['maintenance_schedule']
    q_k = master_result['capacities']
    
    # Solve subproblems in parallel
    subproblem_results = parallel_map(
        solve_subproblem, 
        [(omega, xi_omega, y_k, q_k, d_bar) for omega in scenarios]
    )
    
    # Generate and add Benders cuts (multi-cut, Magnanti-Wong strengthened)
    for omega, result in enumerate(subproblem_results):
        cut = generate_strengthened_cut(result, y_k, q_k, omega)
        master_model.add_cut(cut)
    
    # Check convergence
    ub = compute_upper_bound(subproblem_results, y_k, q_k)
    lb = master_result['objective']
    if (ub - lb) / abs(ub) < tolerance:
        break

5. Case Studies

5.1 Case Study 1: Automotive Coating Facility

5.1.1 System Description

The first case study is based on a representative automotive OEM coating facility producing passenger vehicle body panels. The coating process involves primer application, basecoat, and clearcoat stages, each generating solvent-laden air streams that feed a centralized recovery network. The primary solvents are xylene (aromatic, relatively high value), MEK (ketone, moderate value), and n-butanol (alcohol, lower value). Three recovery units are considered: an activated carbon adsorption system (Unit A), a condenser train (Unit B), and a secondary distillation column for solvent separation and upgrading (Unit C). The planning horizon spans 24 weekly periods, with maintenance windows possible at any period subject to the crew availability constraint of Equation (12) with .

Feed composition uncertainty is substantial: xylene concentration in the exhaust air varies by ±40% around its mean of 2.8 g/m³ depending on the current paint production program. MEK and n-butanol concentrations vary by ±25% and ±30%, respectively. Recovered solvent prices are modeled with weekly autocorrelation of 0.72 (fitted to 36 months of market data), reflecting the persistence typical of petrochemical commodity markets.

5.1.2 Degradation Parameter Estimation

The activated carbon unit (Unit A) exhibits the most significant degradation dynamics. Breakthrough curve data from 18 months of plant operation was used to estimate the degradation rate at nominal throughput (1.0 m³/s exhaust flow) using maximum likelihood estimation with a five-state Markov chain (). The degradation sensitivity coefficients were , , , indicating that the high-boiling-point contaminants in the xylene-rich streams preferentially occupy adsorption sites relevant to all species. Unit B (condenser) shows slower degradation () primarily through fouling film accumulation, while Unit C (distillation) is the most robust () reflecting primarily mechanical wear in the reboiler.

5.1.3 Results

Table 1 summarizes the performance of four decision-making approaches across the 24-week horizon: (a) the deterministic model using mean uncertainty values and ignoring degradation (Baseline); (b) the deterministic model with degradation dynamics but nominal uncertainty (Det+Deg); (c) the two-stage stochastic program without degradation coupling (SP-NoDeg); and (d) the full proposed methodology (SP+Deg). All cost figures are normalized to the Baseline value of 1.000.

The proposed SP+Deg model achieves an 11.9% reduction in expected operating cost relative to the Baseline and, more strikingly, a 22.5% reduction in the CVaR—the tail-cost metric most relevant to operators facing irregular but severe regulatory enforcement actions. The improved CVaR performance reflects the model's ability to schedule proactive maintenance before degradation reaches states where recovery efficiency drops below regulatory thresholds under adverse uncertainty realizations. By contrast, the SP-NoDeg model, while capturing uncertainty in prices and feed compositions, still schedules maintenance suboptimally because it treats equipment efficiency as constant, leading to efficiency surprises in scenarios with high throughput demand.

5.2 Case Study 2: Pharmaceutical API Manufacturing

5.2.1 System Description

The pharmaceutical case study is inspired by a multiproduct API synthesis facility where acetonitrile (MeCN) and dichloromethane (DCM) are the primary process solvents. Both are subject to strict regulatory limits under ICH Q3C guidelines (International Council for Harmonisation, 2021) as Class 2 solvents, meaning recovered solvent must meet purity specifications before reuse in synthesis. The recovery network consists of four units: a distillation column for MeCN recovery (Unit D), a liquid-liquid extraction–distillation hybrid for DCM recovery (Unit E), an activated carbon polishing bed for trace contaminant removal from recovered streams (Unit F), and a thin-film evaporator for concentration of dilute recovery streams (Unit G).

The pharmaceutical context introduces several complications absent from the coating case: (a) purity constraints on recovered solvents are stringent (MeCN ≥ 99.5 mol%, DCM ≥ 99.0 mol%) and are coupled to quality control batch testing, introducing discrete failure modes; (b) the production schedule is irregular, with batch campaigns creating episodic high-solvent-demand periods; (c) regulatory submissions require documentation of solvent recovery compliance, creating long-term record-keeping obligations that penalize retrospective violations more severely than prospective ones.

5.2.2 Uncertainty Sources and Scenario Structure

Feed uncertainty in this context is driven primarily by batch-to-batch variability in reaction yield (affecting solvent consumption) and reaction byproduct generation (affecting recovered solvent impurity profiles). Price uncertainty is less severe than in the coating case—internal transfer pricing for pharmaceutical solvents tends to be more stable—but utility cost uncertainty is significant because the facility is a major electricity consumer. A 250-scenario set is used here, with the scenario tree structured to reflect the batch production calendar: each of the 12 monthly planning periods corresponds to 2–4 batch campaigns with known nominal schedules but uncertain actual yields.

5.2.3 Results and Regulatory Compliance Analysis

The pharmaceutical case highlights a different aspect of the proposed methodology: the interaction between degradation dynamics and purity constraints. As Unit F (the polishing bed) degrades, its ability to remove trace impurities from the recovered MeCN stream declines according to Equation (2), with —a high sensitivity reflecting that trace contaminant removal is the primary function of this unit and is most susceptible to capacity loss. Under the Baseline model (deterministic, no degradation), the optimization allocates high throughput to Unit F in months with intensive batch campaigns, inadvertently accelerating degradation and causing purity failures in months 8–10 when degradation reaches critical levels.

The SP+Deg model, by contrast, predicts this trajectory and schedules a partial regeneration of Unit F at month 6—earlier than the purely economic optimum would suggest—specifically to maintain a safety margin in purity performance during the anticipated high-demand period in months 8–9. This proactive maintenance costs approximately €18,500 in regeneration materials and lost throughput but avoids an estimated €210,000 in batch rejection costs and regulatory penalty exposure (based on the scenario-weighted penalty term), a ratio consistent with real pharmaceutical solvent management economics (Dunn et al., 2004).

Notably, the SP+Deg model achieves 97.6% MeCN purity compliance compared to 91.3% under the Baseline—a substantial improvement with direct implications for batch acceptance rates and ultimately for product yield. The reduction in unplanned shutdowns from 3.2 to 0.6 per year is also economically significant: each unplanned shutdown in this facility incurs estimated costs of €45,000–€80,000 in emergency maintenance labor, replacement materials, and production schedule disruption.

6. Discussion

6.1 Methodological Contributions and Practical Implications

The central methodological contribution of this work is the tight integration of stochastic programming and equipment degradation dynamics within a single optimization model. Previous approaches either optimized operations assuming fixed equipment efficiency (ignoring degradation) or optimized maintenance schedules given fixed operational policies (ignoring operational flexibility). By coupling these decisions through the state evolution equations (Equations 3–4 and 10), the proposed framework reveals a class of solutions—proactive maintenance timed to protect high-value operational windows—that neither approach can discover independently.

This coupling has an important practical consequence for sustainable manufacturing. The SP+Deg model systematically achieves higher recovery rates by keeping equipment closer to its design efficiency envelope; across both case studies, recovered solvent volumes increased by 8–14% relative to the Baseline, directly reducing raw solvent consumption and associated emissions. This aligns with the argument, advanced by Capello et al. (2007) and others, that sophisticated operational optimization of existing recovery infrastructure can yield environmental benefits comparable to capital investment in new separation capacity—at a fraction of the cost.

The risk measure formulation in Equation (5) deserves separate comment. The -weighted CVaR objective allows operators to continuously adjust the risk–performance trade-off as business conditions change. In the pharmaceutical case, where regulatory risk is the dominant concern, a higher (0.7 in our base case) is appropriate. In the coating facility, where margin pressure is the primary driver, a lower (0.3) was preferred by facility management. Sensitivity analysis (not shown in detail here) confirms that the optimal network configuration and maintenance schedule are relatively robust to the choice of over the range 0.2–0.8, with the main impact being on the conservatism of flow routing decisions in high-uncertainty scenarios.

6.2 Scenario Generation and Model Uncertainty

The quality of the stochastic optimization output depends critically on how well the discrete scenario set captures the underlying uncertainty distribution. The moment-matching LHS approach of Høyland and Wallace (2001) performs well for low-dimensional uncertainty (say, up to 10–15 uncertain parameters), but as the number of uncertain parameters grows, moment-matching becomes computationally challenging and the required number of scenarios for adequate representation increases rapidly. In the automotive case, with 9 uncertain parameters (3 compositions × 3 solvents, 3 prices, 1 utility cost), 500 scenarios proved sufficient based on out-of-sample testing. For networks with significantly larger numbers of uncertain streams—common in integrated chemical complexes—scenario reduction techniques (Heitsch & Römisch, 2003) or data-driven scenario generation using generative adversarial networks (Veeramsetty et al., 2022) may be necessary.

The degradation model itself introduces a layer of model uncertainty: the Markov chain structure, the functional form of Equation (2), and the estimated parameters and all represent simplifications and approximations. The linear degradation-efficiency relationship (Equation 2) is a deliberate modeling choice that preserves MILP structure; real systems may exhibit nonlinear relationships, particularly near end-of-life states. Sensitivity analysis on (±30% perturbation) showed that the qualitative recommendations—proactive maintenance, risk-adjusted flow routing—are robust, though the precise timing of maintenance windows shifts by ±1–2 periods.

6.3 Comparison with Alternative Approaches

It is worth comparing the proposed framework against two alternative approaches that are gaining traction in the PSE community. The first is robust optimization (Ben-Tal et al., 2009), which optimizes against the worst case within an uncertainty set rather than over a probability distribution. Robust approaches are computationally lighter and require no distributional information, but they tend to be overly conservative for the types of uncertainty encountered in solvent recovery—feed composition fluctuations are bounded but not adversarially correlated with price movements, so worst-case joint realizations are both extremely rare and excessively penalized. The CVaR-weighted stochastic program of this work provides a more calibrated treatment of tail risk.

The second alternative is data-driven model predictive control (MPC) with a degradation observer, as explored for general chemical processes by Alanqar et al. (2017). MPC is excellent for real-time operational adjustment on short horizons (minutes to hours) but is not naturally suited to the simultaneous optimization of maintenance schedules over weeks-to-months horizons. The rolling horizon implementation described in Section 4.2 can be viewed as a bridge between these paradigms: the stochastic program provides the strategic maintenance schedule and capacity allocation, while a within-period MPC layer (not formalized here) handles fine-grained flow control. Integrating these two layers more tightly—for example, via approximate dynamic programming—is a natural direction for future work.

6.4 Scalability and Computational Considerations

The solve times reported in Table 1—under 50 minutes for the automotive case—are acceptable for weekly replanning cycles but may be limiting for more complex networks. Scalability analysis (not presented in full here) suggests that the L-shaped method's iteration count grows approximately as with network size and linearly with the number of scenarios, consistent with reported behavior for similar decomposition approaches (Birge & Louveaux, 2011). For networks with units, nested decomposition or column generation approaches may be required, or alternatively, scenario reduction to scenarios with compensating out-of-sample robustness assessment.

The implementation was not systematically optimized for parallelism; all 64 available cores were utilized for subproblem solves, but the master program was solved on a single core. Future implementation using parallel master updates (e.g., the asynchronous bundle method) could substantially reduce wall-clock time for larger instances.

7. Conclusion

This article has presented a comprehensive methodology for optimizing solvent recovery networks under uncertainty in feed compositions, market prices, and equipment performance, with explicit coupling of operational decisions to equipment degradation dynamics through a two-stage stochastic programming framework. The key innovations relative to prior work are: (1) the degradation-coupled recovery efficiency model (Equations 2 and 11) that makes equipment health a first-class decision variable; (2) the rolling horizon implementation with sensor-updated degradation state estimates (Equation 15) that enables real-time operational adjustment; and (3) the risk-adjusted CVaR objective (Equations 5–6) that provides a principled approach to balancing expected cost against regulatory and operational tail risk.

Two industrial case studies—an automotive coating facility and a pharmaceutical API manufacturing plant—demonstrate that the proposed approach delivers cost reductions of 12–19% relative to deterministic baselines while substantially improving regulatory compliance rates and reducing unplanned equipment shutdowns. These performance improvements are not achieved through any single mechanism but through the coordinated effect of better-timed maintenance, more flexible flow routing under uncertainty, and explicit protection of high-efficiency operating windows through proactive equipment care.

The implications for sustainable manufacturing are significant. The ability to systematically improve solvent recovery rates through optimization—without capital investment in new equipment—directly reduces raw material consumption and atmospheric emissions. For industries such as pharmaceuticals and coatings that face increasing regulatory scrutiny of VOC emissions and solvent waste streams, the methodology developed here provides a practical pathway to improved environmental performance that is simultaneously economically motivated.

Several directions for future research are evident. First, the integration of this strategic-level stochastic program with a real-time MPC layer, sharing degradation state estimates, would create a complete hierarchical decision support system. Second, the extension to multi-site networks with shared solvent recovery infrastructure—relevant for large chemical parks—introduces network topology decisions and inter-site logistics that substantially expand the decision space. Third, the data-driven estimation of degradation parameters from sparse operational data, using Bayesian inference or transfer learning from similar units, would improve the practical deployability of the approach in settings where extensive historical data is unavailable. These extensions are the subject of ongoing work in our research group.

References