Dynamic River Allocation with GSPy and PuLP Linear Programming (LP)

This example demonstrates how to use the GSPy bridge to dynamically solve a complex, priority-based water allocation problem at every time step of a GoldSim simulation using the external PuLP Linear Programming (LP) library in Python.

Simulating a river network with sequential diversions, return flows, and competing demands governed by priorities can be challenging with traditional explicit rules-based logic. By moving the allocation decision to an external optimization solver, we can achieve optimal water delivery that maximizes overall goal attainment while strictly respecting flow limits and priorities.

River Network Structure

The case study involves a river system with three reaches (N1, N2, N3) and three diversion points (Div A, Div B, Div C).

1. Sequential Flow: The river flows from N1 → N2 → N3 → Sink.
2. Return Flows (R_X): Diversions Div A and Div B have consumption efficiencies (E_A, E_B). The unconsumed portion (Div_X × (1.0 − E_X)) returns downstream:
   • R_A feeds into Reach N2.
   • R_B feeds into Reach N3.
3. Inflows (T_X): Tributaries T2 and T3 add flow at Reaches N2 and N3, respectively.
4. Priorities: Demands must be satisfied in the following order:
   • 1st (Highest): Div B (Weight 3.0)
   • 2nd (Medium): Div A (Weight 2.0)
   • 3rd (Lowest): Div C (Weight 1.0)

When inflows are scarce, the Linear Program simultaneously determines the optimal flows for Div A, Div B, and Div C that satisfy the maximum possible demand while respecting physical flow limits and the priority order.

GSPy and PuLP Implementation

The GoldSim External element utilizes the GSPy bridge to send the current state variables to the Python script, river_optimizer.py. This script then forms the allocation problem and solves it using the PuLP library.

GSPy Interface Definition

The GoldSim External element sends 8 scalar inputs and receives 3 scalar outputs. The order must match the GSPy JSON configuration.

The PuLP Linear Program (LP) Formulation

The problem is formulated as a Weighted Deficit Minimization LP.

1. Decision and Shortage Variables

• Decision Variables: Div_A, Div_B, Div_C (The flows we are solving for, ≥ 0).
• Shortage Variables: S_A, S_B, S_C (The deficit between the diverted flow and the target, ≥ 0).
• Exit Flow: Q_Exit (The flow remaining after Div C, ≥ 0).

2. Objective Function (Minimization)

To enforce the priority order (B > A > C), weights are assigned to the shortage variables. The objective minimizes the weighted shortage, forcing the highest priority demands to be satisfied first.

Minimize (3.0 × S_B) + (2.0 × S_A) + (1.0 × S_C) + (0.01 × Q_Exit)

3. Constraints (Demand & Mass Balance)

The constraints ensure the solution respects both the demand targets and the physical flow limits:

• Demand Constraints: The diverted flow plus any shortage must equal the target:
  • Div_A + S_A = Target_A
  • Div_B + S_B = Target_B
  • Div_C + S_C = Target_C
• Reach Limits (N1): Div_A cannot exceed T1.
  • Div_A ≤ T1
• Reach Limits (N2): Flow available at N2 is the remainder from N1 plus tributary T2 and return flow R_A (R_A = Div_A × (1.0 − E_A)).
  • Div_B ≤ (T1 − Div_A) + T2 + R_A
• Reach Limits (N3): Flow available at N3 is the remainder from N2 plus tributary T3 and return flow R_B (R_B = Div_B × (1.0 − E_B)).
  • Div_C ≤ (N2 Flow − Div_B) + T3 + R_B
• Exit Flow Balance: Defines the final flow to the sink.
  • Q_{in\_N3} − Div_C = Q_Exit

Contact:

Jason Lillywhite (Lillywhite Water Solutions)

Conclusion and Download

This integration between GoldSim, GSPy, and PuLP allows the GoldSim model to ask, "What is the best allocation?" at every time step, moving from a rigid rules-based simulation to a true dynamic optimization.

This approach can be adapted to solve various resource allocation problems within the robust probabilistic framework of GoldSim.

Download the complete example package below.

Setup Note: You must have a 64-bit Python installation (3.8+ recommended) and the NumPy and PuLP packages installed. The downloaded DLL and JSON file must have the same names and be in the same directory to function correctly. Please refer to the GSPy documentation for more instructions on configuration and troubleshooting.