Lagrangian Decomposition and Parallel Optimization for Integrated Pricing and Replenishment

1. Business Problem

Retail and distribution planning decisions are strongly connected. A price change affects demand; demand affects inventory depletion; inventory requirements affect purchase quantities and timing; and purchases consume supplier and distribution-center capacity. Solving these decisions independently can create plans that appear attractive locally but are infeasible or expensive when evaluated as one operating plan.

The model supports tactical planning decisions such as:

• which price option to use for each SKU, distribution center, and planning period;

• how much of the price-dependent demand to serve;

• how much to purchase from each eligible supplier and when to place each order;

• which distribution center should receive each order, with arrival timing determined by the supplier-to-DC lead time;

• how much ending inventory to hold by SKU, distribution center, and period;

• how to coordinate these decisions while satisfying supplier-availability, receiving-capacity, and storage-capacity constraints. 

The formulation provides a tactical planning tool for coordinating pricing, sales, purchasing, inbound flows, and inventory across suppliers, distribution centers, and planning periods. It can also support scenario analysis involving tighter capacity, supplier disruptions, alternative lead times, different price ladders, and changes in price-dependent demand. 

2. Decision Scope and Solution Framework

The planning system connects five decisions across products, suppliers, distribution centers, and periods: price selection, sales, purchasing, inbound timing, and inventory. Price-dependent demand determines potential sales, supplier eligibility and availability constrain replenishment, lead times determine arrival periods, and DC receiving and storage capacities couple all SKU plans. The solution framework below compares the integrated exact MILP with the SKU-decomposed Lagrangian approach, including parallel subproblem execution, feasibility repair, multiplier updates, and stopping evaluation.

3. Mathematical Model

The integrated planning problem is formulated as a mixed-integer linear program. Binary variables determine the selected price option for each product, distribution center, and planning period. Continuous variables represent sales, supplier purchases, and ending inventory. Inbound arrivals are derived from purchase decisions, supplier-to-DC lead times, and initial pipeline quantities.

3.1 Sets and Indices

3.2 Parameters

3.3 Decision variables and Derived Arrivals 

3.4 Objective Function

The objective maximizes total operating profit, defined as sales revenue minus promotion, purchasing, transportation, and inventory-holding costs:

3.5 Constraints

4. Motivation for Lagrangian Decomposition

The exact MILP grows quickly because price-selection binaries and flow variables are indexed by product, location, period, and, for purchases, supplier. Most constraints are product-specific, but DC receiving and storage capacities couple all SKUs. These shared constraints prevent the SKUs from being solved independently and therefore become the main target for decomposition.

The computational strategy separates the model by SKU while retaining a mechanism that prices shared-capacity consumption. This creates many smaller subproblems that can be solved sequentially or in parallel.

5. Lagrangian Relaxation

5.1 Relaxed Constraints and Penalized Objective 

5.2 Sequential and Parallel Subproblem Execution 

The relaxed model contains one independent subproblem for each SKU. In sequential mode, these subproblems are solved one after another. In parallel mode, the SKU set is distributed across multiple worker processes, and each worker solves its assigned subproblems independently using the same Lagrangian multipliers.

After all workers finish, their SKU-level solutions, objective values, and capacity usage are collected and aggregated before violation measurement, repair, bound evaluation, and multiplier updates. Therefore, parallel execution changes only how the subproblems are solved; it does not change the mathematical formulation, repair procedure, bounds, or stopping criteria.

The implementation uses six worker processes with one solver thread per worker. This avoids excessive nested parallelism while allowing many independent SKU subproblems to be solved concurrently.

5.3 Per-iteration Solution Flow 

Each Lagrangian iteration follows the same sequence:

• solve one independent pricing-and-replenishment subproblem for each SKU using the current multipliers;

• aggregate SKU decisions into a complete relaxed portfolio;

• calculate total DC receiving and storage usage and identify violations;

• evaluate the relaxed objective and update the best dual value;

• if the aggregated solution violates shared capacity, invoke the repair heuristic to obtain an original-model feasible solution;

• evaluate the repaired solution as a primal incumbent;

• update the multipliers based on capacity violations and continue until the gap, iteration, no-improvement, or time criterion is met.

5.3 Multiplier Update

The initial step is small to avoid destabilizing near-feasible solutions. The step may grow within a configured limit when violations persist and can be reduced when the method approaches feasibility or when progress stalls. The implementation records the current dual value, best dual value, current and best primal values, maximum violation, feasibility state, step size, iteration time, and no-improvement count.

5.4 Bounds and Stopping

For a maximization problem, the relaxed solution provides an upper bound when subproblem bounds are valid, while any repaired solution that satisfies the original constraints provides a lower bound. The relative bound gap is monitored as:

Gap = (Best Dual − Best Primal) / max{|Best Primal|, ε}

The procedure stops when the acceptable gap is reached, the iteration or time limit is reached, or no meaningful improvement is observed for the configured number of iterations. Sequential and parallel modes execute the same mathematical method; only the subproblem execution strategy changes.

6. Capacity Repair Heuristic

The relaxed portfolio can violate shared DC capacities because those constraints are not enforced inside individual SKU subproblems. The repair method converts that portfolio into a feasible operating plan while attempting to preserve as much objective value as possible.

6.1 Search Structure

The repair is a bounded, deterministic local search. It ranks active receiving and storage violations by excess amount, identifies the SKUs contributing most to the selected violation, generates a limited set of candidate moves, scores them, applies the best admissible move, and repeats until capacity is feasible or a configured limit is reached.

Candidate phases are evaluated in an operationally conservative order:

· period shift: move inbound quantity to a nearby period at the same DC;

· DC reallocation: redirect quantity to another DC when timing and capacity permit;

· supplier substitution: move quantity to another eligible supplier and arrival pattern;

· purchase reduction: reduce inbound quantity when relocation options are unavailable;

· price change: use a demand-increasing price option as a final storage-relief mechanism when excess inventory is the remaining issue.

6.2 Candidate Evaluation and Optimization

To avoid repeatedly copying and recalculating the full solution, each candidate uses a sparse copy-on-write overlay. Unchanged decisions are read from the current solution, while only modified values are stored in the candidate patch. Inventory is recomputed only from the earliest affected period until the trajectory rejoins the unchanged base solution.

Capacity and supplier usage are maintained incrementally. Candidate scoring favors moves that reduce more joint capacity violation with less objective loss and less lost revenue. A simplified representation is:

Score = (Violation Reduction − Newly Created Violation) / (Objective Loss + Lost Revenue + ε)

Here, ε is a small positive constant used to prevent division by zero when the objective loss and lost revenue are both zero. Only moves that produce a positive net reduction in total receiving and storage violation are eligible. Before declaring success, the cached capacity state is reconciled against a full recalculation, and the final plan is checked against every original-model constraint.

7. Implementation Architecture

The implementation follows a modular optimization-software architecture that separates data generation, mathematical-model construction, solver execution, Lagrangian coordination, feasibility repair, validation, visualization, and experiment reporting. The exact MILP and SKU-decomposed approaches use the same generated problem instance and economic definitions, ensuring that differences in objective value, feasibility, and runtime arise from the solution method rather than inconsistent data or model logic.

The architecture supports four execution paths:

• exact MILP with Gurobi;

• exact MILP with HiGHS;

• sequential SKU-decomposed Lagrangian relaxation;

• parallel SKU-decomposed Lagrangian relaxation.

Both solver families implement the same underlying formulation. This provides an independent cross-solver validation mechanism while preserving solver-specific model construction and execution logic.

7.1 Shared Data and Configuration Layer

pricing_replenishment_config.py provides the central experiment configuration, including problem dimensions, random seed, solver controls, Lagrangian settings, parallel-worker configuration, stopping criteria, repair limits, and output options.

pricing_replenishment_data.py generates the complete deterministic planning instance, including:

• products, suppliers, distribution centers, periods, and price options;

• eligible product-supplier relationships;

• price-dependent demand;

• purchase, transportation, promotion, and holding costs;

• supplier availability;

• supplier-to-DC lead times;

• initial inventory and pipeline arrivals;

• DC receiving and storage capacities;

• minimum-inventory requirements.

A fixed random seed allows the same instance to be regenerated across exact, sequential, and parallel experiments. This shared data layer is critical for valid comparisons because every solution method receives the same products, costs, capacities, demand values, and network relationships.

pricing_replenishment_common.py contains reusable structures and utility logic shared across the solution workflows. pricing_replenishment_solution.py provides standardized solution representations so exact and Lagrangian methods can return comparable decisions and performance measures.

7.2 Mathematical Model Layer

pricing_replenishment_model.py contains the Gurobi formulation used by both the integrated exact model and the product-restricted Lagrangian subproblems. The same model-building function can construct:

• the original integrated MILP with shared DC-capacity constraints; or

• a restricted SKU subproblem in which those constraints are removed and capacity usage is priced through Lagrangian multipliers.

This design avoids maintaining separate mathematical formulations for the exact and decomposed approaches. SKU-local constraints including price selection, demand linkage, supplier availability, lead-time-adjusted arrivals, inventory balance, and variable bounds remain identical in both cases.

The HiGHS implementation is maintained in the highs_solve package. Its model and execution components mirror the Gurobi structure through:

• pricing_replenishment_highs_model.py;

• pricing_replenishment_highs_exact.py;

• pricing_replenishment_highs_lagrangian.py;

• pricing_replenishment_highs_lagrangian_parallel.py.

Keeping the HiGHS implementation separate from the Gurobi API prevents solver-specific syntax from contaminating the shared application layer while retaining mathematical consistency between the two formulations.

7.3 Exact and Decomposed Solution Engines

pricing_replenishment_exact.py executes the integrated Gurobi MILP. The corresponding HiGHS engine solves the same planning problem through the HiGHS model interface. Both exact workflows return standardized information such as objective value, best bound, relative gap, solver status, runtime, and solution availability.

pricing_replenishment_lagrangian.py implements the sequential SKU-decomposed procedure. For each iteration, it:

1.     solves one independent subproblem for every SKU;

2.     aggregates the SKU-level decisions;

3.     calculates shared receiving and storage usage;

4.     evaluates the dual bound;

5.     invokes the repair heuristic when the relaxed portfolio is infeasible;

6.     evaluates the resulting primal solution;

7.     updates the capacity multipliers;

8.     checks the stopping criteria.

pricing_replenishment_lagrangian_parallel.py implements the same mathematical procedure using multiple worker processes. SKU subproblems are distributed across workers and solved concurrently using identical multipliers. Their solutions, bounds, and capacity contributions are then returned to the parent process and aggregated before repair, multiplier updates, and stopping evaluation.

The sequential and parallel engines therefore differ only in subproblem execution. They use the same formulation, multiplier logic, repair method, bound calculations, and termination rules. This makes runtime comparisons meaningful and allows parallel correctness to be verified by comparing complete iteration histories and final solutions.

7.4 Capacity-repair Layer

pricing_replenishment_repair.py contains the unified feasibility-repair procedure used by both sequential and parallel Lagrangian workflows. The repair logic is separated from the decomposition engine so that it can be tested, optimized, and extended independently.

The implementation uses a bounded deterministic local search with several candidate-move families:

• period shifting;

• distribution-center reallocation;

• supplier substitution;

• purchase reduction;

• price-option adjustment.

To support large instances, candidate evaluation uses sparse copy-on-write modifications, incremental capacity tracking, cached supplier usage, and partial inventory recomputation from the earliest affected period. These mechanisms avoid repeatedly copying and recalculating the full portfolio for every candidate move.

The repair component also records detailed diagnostics, including the number and type of accepted moves, objective loss, lost sales quantity, remaining violation, stopping reason, and candidate counts by repair phase. This makes the heuristic behavior observable rather than treating it as a black-box feasibility adjustment.

7.5 Orchestration and Execution Flow

pricing_replenishment_pipeline.py coordinates the complete experiment lifecycle:

1.     validate the configuration;

2.     create a unique experiment identity;

3.     generate the deterministic planning instance;

4.     invoke the selected solver and solution method;

5.     validate the returned solution against the original model;

6.     generate diagnostic plots;

7.     append standardized records to shared output files.

pricing_replenishment_run.py provides a focused execution entry point, while pricing_replenishment_progress.py standardizes timestamps, stage reporting, elapsed-time measurements, and iteration progress.

The pipeline records both solver-runtime sums and actual algorithm wall time. This distinction is particularly important for parallel Lagrangian runs because the sum of individual subproblem runtimes may exceed elapsed wall time when several subproblems are solved concurrently.

Each run receives a unique identifier containing the solver, method, execution mode, problem dimensions, and random seed. This identity is used consistently across logs, plots, and CSV records, allowing all generated artifacts to be traced back to their exact configuration.

7.6 Validation and Objective Consistency

pricing_replenishment_feasibility.py independently checks the returned solution against the original integrated constraints. The validation covers:

• DC storage capacity;

• DC receiving capacity;

• supplier availability;

• inventory balance;

• demand and price linkage;

• nonnegative quantities.

This validation is applied to exact solutions and repaired Lagrangian solutions. A decomposed solution is therefore not accepted merely because the repair procedure reports success; it must also pass a complete original-model feasibility check.

The implementation additionally recalculates the original economic objective from the returned decisions and compares it with the aggregated subproblem objective. This detects inconsistencies in decision aggregation, Lagrangian accounting, repair updates, and solver-result extraction.

Cross-solver and cross-execution validation provides additional safeguards:

• Gurobi and HiGHS exact solutions are compared on small instances;

• sequential and parallel Lagrangian modes are expected to reproduce the same objective, bound, repair behavior, and stopping condition;

• exact and decomposed methods are evaluated using the same data and original objective definition.

7.7 Reporting and Experiment Output

pricing_replenishment_export.py writes standardized experiment records to shared CSV files. Depending on the selected workflow, outputs include:

• run-level objective, bound, gap, runtime, and feasibility;

• iteration-level dual, primal, violation, step-size, and timing history;

• selected prices, sales, and inventory by product, DC, and period;

• supplier purchases and lead-time-adjusted arrivals;

• DC receiving- and storage-capacity utilization;

• supplier-availability utilization;

• repair diagnostics and termination information.

pricing_replenishment_plot.py generates diagnostic figures for:

• Lagrangian primal and dual progression;

• shared-capacity violations;

• multiplier step sizes;

• DC storage utilization;

• DC receiving utilization.

Generated outputs are stored separately from the source files in the outputs directory. Transient Python artifacts and generated results are excluded from version control through .gitignore, keeping the repository focused on source code, configuration, documentation, and reproducible execution logic.

7.8 Software Architecture and Design Strengths

The implementation demonstrates several software-engineering strengths that support maintainability, validation, and reproducible experimentation.

Single source of model logic.
The integrated Gurobi model and the Gurobi SKU subproblems are generated from the same model-building component, reducing duplicated logic and formulation drift.

Separation of responsibilities.
Data generation, model construction, solver execution, Lagrangian coordination, repair, feasibility validation, plotting, and export are implemented as distinct components with clearly defined responsibilities.

Solver-specific isolation.
Gurobi and HiGHS use separate solver-specific model layers while preserving the same mathematical formulation. This prevents solver APIs from being mixed with the shared application logic.

Reusable solution structures.
Exact, sequential Lagrangian, and parallel Lagrangian workflows return standardized solution objects and common performance measures, simplifying validation, reporting, and method comparison.

Independent validation layer.
Returned solutions are checked outside the optimization and repair routines against the original model constraints. This separates solution generation from solution verification.

Reproducible experiment management.
Centralized configuration, fixed random seeds, unique run identifiers, consistent timing, and structured output files allow experiments to be repeated and traced to their exact settings.

Modular parallel execution.
Parallel worker management is separated from the mathematical subproblem definition. The same SKU subproblem logic can therefore be used in sequential or parallel execution.

Extensibility.
New solver backends, constraints, solution methods, diagnostics, or reporting outputs can be added without restructuring the complete application.

Table 1 summarizes the main repository files and their responsibilities.

7.9 Repository and Reproduciblity

The complete implementation is maintained in the following public repository:

GitHub repository:
https://github.com/asefi-h/pricing-replenishment-optimization

The repository contains the Gurobi and HiGHS implementations, exact and SKU-decomposed solution engines, configuration files, dependency specification, execution instructions, and project documentation.

Reproducibility is supported through:

• centralized experiment configuration;

• fixed random seeds for deterministic instance generation;

• shared data and objective definitions across all solution methods;

• standardized run identifiers and CSV outputs;

• documented solver settings and execution commands;

• version-controlled source code and commit history;

• a GitHub Actions continuous-integration workflow;

• an automatically updated CI status badge in the repository README.

The CI workflow runs for pushes and pull requests targeting the main branch. It prepares a clean Ubuntu environment, configures Python 3.12, installs the dependencies listed in requirements.txt, verifies that the main pipeline can be imported, and compiles the Python source files.

The CI process provides a lightweight check of environment setup, dependency installation, imports, and source-code syntax. Mathematical validation, solution feasibility, objective consistency, and scalability are evaluated separately through the computational experiments and independent validation routines described in this study.

8. Computational Results and Scalablity

The computational study compares the integrated exact MILP with sequential and parallel SKU-decomposed Lagrangian relaxation using HiGHS. Runtime is reported as final algorithm wall time.

For the smallest instances, parallel execution is slower because process-management overhead exceeds the benefit of concurrent subproblem solution. The crossover occurs at 500 SKUs, where parallel Lagrangian relaxation is 2.20 times faster than sequential execution and 5.49 times faster than the exact MILP.

The advantage becomes substantial at 1,000 SKUs. Parallel Lagrangian relaxation completes in 19.86 seconds, compared with 56.33 seconds sequentially and 643.88 seconds for the exact MILP. This corresponds to a 2.84-fold speedup over sequential decomposition and a 32.41-fold speedup over the integrated exact solve.

At 2,000 and 5,000 SKUs, the exact MILP reaches the time limit without producing an incumbent. In contrast, both Lagrangian approaches return solutions that satisfy the original model constraints. Parallel execution completes these cases in 62.72 and 131.19 seconds, respectively.

The method continues to scale to 10,000 and 20,000 SKUs. Parallel execution completes the 10,000-SKU case in 296.46 seconds and the 20,000-SKU case in 796.60 seconds, providing speedups of 3.02 and 2.51 relative to sequential execution. Across all reported Lagrangian cases, the retained solutions are feasible in the original integrated formulation, and the final relative gaps remain below 0.1%.

Overall, the results show a clear computational transition: the integrated MILP is competitive only for small instances, while parallel SKU decomposition provides the strongest runtime performance and remains effective at problem sizes for which the exact formulation does not produce a usable solution within the recorded solve horizon.

Lagrangian Bound Convergence: The representative 300-SKU run illustrates the progression of current and retained primal and dual values. The best feasible lower bound improves materially at iteration 8, while the best Lagrangian upper bound continues to decrease. By iteration 12, the retained bounds are nearly coincident and the acceptable-gap stopping criterion is satisfied.