Large-Scale Linear Programming for Urban Water Network Optimization Using Interior Point Methods

The objective of this study is to optimize the allocation of clean water from multiple water treatment plants to various demand zones across Medan city. The problem is modeled as a large-scale linear program that considers production capacities, consumer demands, operational costs, hydraulic pressure constraints, and transmission losses. The distribution network includes 6 major water treatment plants (IPA) and 21 distribution zones served by PDAM Tirtanadi. Key operational challenges include excessive NRW, inefficient routing, and pressure disparities across different elevation zones.

2.1 Problem description

PDAM Tirtanadi is the regional water utility responsible for supplying clean water to over 2.4 million residents in the city of Medan, Indonesia. The distribution system spans 21 service zones and is supported by multiple IPAs, including major facilities in Sunggal, Deli Tua, Martubung, Denai, and others. Each plant has a fixed daily production capacity and serves a subset of demand zones, often located at varying elevations. The network is challenged by operational inefficiencies, including NRW losses of up to 26.5%, underutilized capacity in certain plants, overburdened transmission routes, and significant pressure drops in elevated zones. Additionally, infrastructure age and physical topography contribute to unequal service reliability across zones.

Recent advancements in infrastructure modeling increasingly rely on hybrid optimization and intelligent systems to address uncertainty, sustainability, and performance trade-offs. Roushan et al. [10] utilized neutrosophic programming with blockchain-based smart contracts to optimize disaster relief supply chains under uncertainty, offering a transparent and adaptive framework. Salamian et al. [11] introduced a robust multi-product bioenergy supply chain model, addressing uncertainty using fuzzy torrefaction parameters. Al Momin et al. [12] applied fuzzy AHP and geospatial modeling to optimally site electric vehicle charging stations, exemplifying hybrid spatial decision-making. Dorahaki et al. [13] proposed a robust optimization model to manage energy flexibility and environmental sustainability in local energy communities. Albahri et al. [14] provided a taxonomy of trustworthy AI applications in natural disaster response, integrating fuzzy logic, explainable AI, and multi-criteria decision-making.

Hafeznia and Stojadinović [15] developed the ResQ-IOS framework, using iterative simulation-optimization for interdependent infrastructure resilience. Chen et al. [16] systematically reviewed multi-objective models in energy system planning, linking long-term sustainability with social and economic criteria. Abualigah et al. [17] examined metaheuristics for sustainable supply chains, highlighting optimization’s role in balancing environmental impact and efficiency. Kakamoukas et al. [18] proposed a fuzzy-logic-driven UAV routing protocol for smart agriculture, enhancing mission efficiency in uncertain environments. Lastly, Zangato et al. [19] introduced an RL-based energy management system integrating clustering and constrained learning to improve adaptability in smart buildings. These studies demonstrate the power of hybrid mathematical tools to manage uncertainty across sectors. Building on this foundation, our work contributes a fuzzy-stochastic-robust LP framework to urban water networks, bridging mathematical theory and operational impact [20, 21].

From an economic perspective, each treatment plant has its own cost structure based on water treatment methods (conventional or membrane-based), energy consumption for pumping, and pipeline maintenance. These costs vary from Rp 2.8 to Rp 3.6 per cubic meter. Transporting water over longer distances or to higher elevation zones incurs additional energy losses, which are typically modeled using head-loss equations. Because operational costs and hydraulic losses are location-specific and interdependent, suboptimal routing can result in disproportionately high costs or service disruptions in peripheral zones.

This study formulates the problem as a large-scale linear optimization model aimed at minimizing the total operational cost while ensuring full demand coverage, respecting supply limitations, and maintaining minimum hydraulic pressure requirements throughout the network [22, 23]. The resulting model reflects the core trade-offs faced by PDAM Tirtanadi: how to allocate limited resources efficiently across a geographically diverse urban service area under both physical and economic constraints.

2.2 Mathematical model formulation

This section formulates the urban water distribution optimization problem as a large-scale LP model [24, 25]. The objective is to minimize the total operational cost of distributing clean water from multiple IPAs to various demand zones, subject to capacity, demand, and hydraulic constraints. To represent the system mathematically, we define a set of decision variables that determine the optimal water flow from each source to each destination. Parameters include water production capacities, zonal demand volumes, cost coefficients, and elevation-based pressure requirements. The constraints in the model ensure that the total water allocated from each IPA does not exceed its capacity, that every demand zone receives at least its required volume, and that minimum hydraulic pressure levels are maintained. Additionally, the model incorporates head loss across pipelines using a linearized approximation of the Hazen-Williams equation, allowing the optimization to remain in LP form. This formulation results in a sparse, high-dimensional optimization problem well-suited for interior point solution techniques.

Objective function

Minimize total operational cost

$\min Z=\sum\limits_{i\in I}^{{}}{\sum\limits_{j\in J}{{{C}_{ij}}}}.{{x}_{ij}}$                       (1)

subject to,

1. Supply capacity constraint (per plant):

$\sum\limits_{j\in J}{{{x}_{ij}}\le {{S}_{i}},\begin{matrix}   {} & {}  \\\end{matrix}\forall i\in I}$                             (2)

2. Demand fulfillment constraint (per zone):

$\sum\limits_{i\in I}{{{x}_{ij}}\ge {{D}_{j}}\begin{matrix}   {} & {}  \\ \end{matrix}\forall j\in J}$                (3)

3. Hydraulic pressure constraint (linearized):

${{H}_{i}}-{{K}_{ij}}.{{x}_{ij}}\ge {{h}_{j}}\begin{matrix}   {} & {}  \\ \end{matrix}\forall i\in I,j\in J$                    (4)

4. Non-negativity constraint:

${{x}_{ij}}\ge 0\begin{matrix}   {} & {}  \\ \end{matrix}\forall i\in I,j\in J$                      (5)

The proposed mathematical model is formulated as a large-scale linear programming problem designed to optimize the operational efficiency of an urban water distribution network. The objective function minimizes the total cost of transporting water from multiple treatment plants to various demand zones, where the unit cost ${{C}_{ij}}$ includes treatment, pumping, and maintenance expenses. The model includes four main constraint groups:

(1) The supply capacity constraint ensures that the total volume of water distributed from each plant does not exceed its production capacity ${{S}_{i}}.$

(2) The demand fulfillment constraint guarantees that each zone receives at least its required daily demand ${{D}_{j}}.$

(3) The hydraulic pressure constraint, represented in linearized form as ${{H}_{i}}-{{K}_{ij}}.{{x}_{ij}}\ge {{h}_{j}}$, ensures sufficient pressure at the consumer end to maintain service quality.

(4) The non-negativity constraint enforces that all water flow variables remain physically meaningful. Together, these components produce a high-dimensional sparse LP model that captures the critical engineering, economic, and hydraulic aspects of water distribution planning, while remaining solvable using efficient interior point algorithms.

To formalize the optimization problem of urban water distribution, a mathematical model was constructed using linear programming principles. The model integrates key system components, including supply sources, demand zones, cost structures, and hydraulic constraints, into a unified decision framework. Each element of the system is translated into mathematical terms sets, parameters, variables, objective function, and constraints to ensure both operational feasibility and cost efficiency. The overall structure of this mathematical approach is illustrated in Figure 1, which outlines the core components of the model and their interrelationships within the optimization process.

image008.png

Figure 1. Mathematical approach

As shown in Figure 1, the model formulation begins by defining the system structure through sets and indices that distinguish water treatment plants and demand zones. Parameters such as production capacities, daily demands, unit costs, and elevation-related pressure requirements are embedded to capture the physical and economic dynamics of the distribution network. The decision variable represents the volume of water transported from each plant to each zone, while the objective function seeks to minimize total operational costs. The constraints ensure that supply limits are not exceeded, demand requirements are satisfied, minimum hydraulic head is maintained, and all variable values remain non-negative. This structured representation enables the model to be solved efficiently using interior point algorithms, which are well-suited for handling the large number of decision variables and sparse constraint matrices typical of real-world urban utility systems.

To model the distribution system accurately, it is essential to understand the structural characteristics and operational capacity of the major water treatment plants managed by PDAM Tirtanadi. Each plant differs in its daily production capacity, the number of zones it serves, the average elevation of those zones, and its associated production costs. These variations significantly impact the optimization model, as they influence both the feasibility of meeting demand and the cost-effectiveness of supply routes. The production cost values shown in Table 1 were sourced from PDAM Tirtanadi’s internal operational records and engineering reports covering the 2023–2024 fiscal period. These figures include unit energy consumption, chemical treatment expenses, and staff overhead, as confirmed by an internal energy audit conducted in June 2023.

Table 1. Capacities and zones

Estimated unit production costs (Rp/m³) by treatment plant. Data sourced from PDAM Tirtanadi operational reports and energy audits (2023–2024) (Table 1).

As shown in Table 1, each water treatment plant differs in its production capacity, zone coverage, and operational costs. Notably, plants with higher capacities such as Sunggal and Deli Tua exhibit lower unit costs per cubic meter due to economies of scale. Conversely, smaller plants such as Martubung and Limau Manis face higher per-unit costs, which may be attributed to shorter pipeline networks or higher energy requirements for localized pumping.

2.3 Linearization techniques

In real-world water distribution systems, head loss due to friction in pipelines is a nonlinear function typically governed by empirical formulas such as the Hazen-Williams or Darcy-Weisbach equations. The Hazen-Williams formula defines head loss as:

${{h}_{L}}=10.67\times \frac{L}{{{C}^{1.85}}.{{d}^{4.87}}}\times {{Q}^{1.85}}$                   (6)

where, L is the pipe length (m), C is the Hazen-Williams roughness coefficient, d is the pipe diameter (m), Q is the discharge (flow rate in m³/s). However, this model is nonlinear and non-convex, making it incompatible with linear programming (LP) solvers. To preserve computational efficiency and maintain compatibility with interior point methods, the nonlinear head-loss term is approximated linearly as,

${{h}_{L}}\approx {{K}_{ij}}.{{x}_{ij}}$

where, ${{K}_{ij}}$ is a linearized head loss coefficient between water treatment plant i and demand zone j, empirically estimated based on pipe length, elevation, material, and average operating conditions. This simplification allows the pressure constraint:

${{H}_{i}}-{{K}_{ij}}.{{x}_{ij}}\ge {{h}_{j}}$

This simplification maintains model solvability while reflecting realistic behavior. To validate this linearization, we compared it against the nonlinear Hazen-Williams formulation using operational data from PDAM Tirtanadi. Over 120 pipeline segments, the linear model yielded an RMSE of 0.73 meters and a MAPE of 5.8%, confirming its appropriateness for strategic-level planning. The resulting LP problem is formulated as:

$\min {{\mathbf{c}}^{T}}\mathbf{x}\begin{matrix}   {} & \text{subject to}  \\\end{matrix}\mathbf{Ax}\le \mathbf{b},\begin{matrix}   {} & \mathbf{x}\ge 0  \\ \end{matrix}$        (7)

where, x is the vector of flow decisions ${{x}_{ij}}$, c is the cost vector composed of ${{C}_{ij}}$, A is the constraint matrix representing supply, demand, and pressure limits, and b is the vector of upper/lower bounds (capacities, demands, minimum pressures). This transformation ensures the problem is solvable using primal-dual interior point methods, which are particularly effective for large-scale sparse LPs such as those encountered in water distribution optimization. Although the actual behavior of urban water networks is governed by unsteady hydraulics and spatial variability (e.g., changes in pipe diameter and valve dynamics), this model is designed for long-term strategic planning, where the assumption of steady-state flow and static geometry is considered appropriate. Transient conditions and local pipe heterogeneity are typically handled in real-time supervisory control systems, not in high-level optimization. The use of averaged diameter classes and nominal operating flow conditions allows the model to remain computationally tractable while still producing meaningful and implementable policy recommendations. Similar simplifications have been adopted in urban infrastructure planning models using linearized hydraulics.

2.4 Linearization techniques

To solve the formulated LP problem:

$\min {{\mathbf{c}}^{T}}\mathbf{x} \text{subject to }\mathbf{Ax=b} \mathbf{x}\ge 0$                (8)

We employ the primal-dual path-following interior point method, which iteratively solves a sequence of perturbed KKT (Karush-Kuhn-Tucker) systems while maintaining the solution in the interior of the feasible region. The dual problem of the LP is given by:

$  \underset{y,z}{\mathop{\max }}\,{{\mathbf{b}}^{T}}\mathbf{y} \text{subject to }{{\mathbf{A}}^{T}}\mathbf{y}+\mathbf{z=c} \mathbf{z}\ge 0$               (9)

The central path condition adds a barrier term:

$ XZ=\mu \mathbf{e}$                     (10)

where, $X=\text{diag}\left( x \right)$, $Z=\text{diag}\left( z \right)$, and μ is the barrier parameter (reduced each iteration).

To find a search direction $\left( \Delta x,\Delta y,\Delta z \right)$, the algorithm solves the Newton system derived from the perturbed KKT conditions,

where, ${{r}_{p}}=b-Ax$ is the primal residual, ${{r}_{d}}=c-{{A}^{T}}y-z$ is the dual residual, ${{r}_{c}}=\mu e-XZ$ is the centering residual (complementarity gap). After solving for the direction, the algorithm updates the variables with step size $\alpha \in 0,1$.

$ \begin{matrix}   {{\mathbf{x}}^{\left( k+1 \right)}} & = & {{\mathbf{x}}^{\left( k \right)}}+\alpha \Delta \mathbf{x}  \\   {{\mathbf{y}}^{\left( k+1 \right)}} & = & {{\mathbf{y}}^{\left( k \right)}}+\alpha \Delta \mathbf{y}  \\   {{\mathbf{z}}^{\left( k+1 \right)}} & = & {{\mathbf{z}}^{\left( k \right)}}+\alpha \Delta \mathbf{z}  \\ \end{matrix} $              (12)

The parameter μ is updated as,

$\mu =\frac{{{\mathbf{x}}^{T}}\mathbf{z}}{n}$                     (13)

And the algorithm proceeds until $\mu <\varepsilon $, where, ε is the convergence tolerance.

2.4 Algorithm implementation in MATLAB

To operationalize the mathematical model, a primal-dual interior point algorithm is employed due to its superior performance in solving large-scale linear programming problems with sparse constraint matrices. Unlike the simplex method, which traverses the edges of the feasible polyhedron, the interior point method follows a trajectory through the interior, guided by the central path defined by a logarithmic barrier function. This approach allows for faster convergence and better numerical stability, particularly in high-dimensional systems such as urban water distribution networks. The algorithm is implemented in MATLAB and structured into three key phases: initialization, iterative update steps based on Newton's method, and solution extraction. The step-by-step procedure is outlined in the interior point algorithm.

The interior point algorithm provides a numerically efficient mechanism to solve the high-dimensional optimization problem underlying the urban water distribution system. By maintaining feasibility at every iteration and simultaneously updating both primal and dual variables, the method guarantees rapid convergence to the global optimum. The use of sparse matrix techniques and Newton-based updates ensures that the algorithm remains scalable even when applied to real-world systems involving dozens of supply sources and distribution zones. The output includes not only the optimal flow values that minimize cost, but also insightful economic indicators such as shadow prices and marginal costs, which can inform strategic decisions in network expansion, tariff setting, and infrastructure prioritization. This algorithmic framework forms the computational core of the water distribution optimization model presented in this study.

To capture the spatial and operational dynamics of the urban water distribution network, a connection matrix was developed to map the relationships between water treatment plants and their corresponding service zones. This matrix reflects actual infrastructure linkages, based on pipeline connectivity, elevation alignment, and operational routing from PDAM Tirtanadi. Each entry indicates whether a treatment plant directly serves a particular service region. In addition to network structure, the table includes critical demand-side attributes such as total population, number of customers, and average water supply (in liters per second) for each zone. These values provide a comprehensive basis for formulating demand constraints and validating flow allocations in the optimization model. The detailed plant-to-zone connections and demand characteristics are summarized in Table 2.

Table 2. Plant-to-zone connection matrix

As shown in Table 2, the distribution of service coverage across plants is heterogeneous. High-capacity plants such as Sunggal and Deli Tua are linked to multiple service regions, including densely populated areas such as Medan Helvetia and Medan Johor. On the other hand, smaller plants like Limau Manis, Martubung, and Pulo Brayan provide localized support to peripheral zones or elevated regions that require additional pumping effort. The variation in population, customer base, and water supply across zones illustrates the underlying complexity of demand allocation. For instance, Medan Helvetia receives 325 L/s with a population base exceeding 150,000, while Medan Petisah is supplied with just 170 L/s for a smaller yet strategically critical zone. These disparities must be carefully addressed within the optimization framework to ensure equitable, cost-effective, and hydraulically feasible water distribution throughout the city. The connection matrix also forms the structural basis for the constraint system used in the linear programming model.