Surrogate “Level-Based” Lagrangian Relaxation for mixed-integer linear programming

Mixed-Integer Linear Programming (MILP) plays an important role across a range of scientific disciplines and within areas of strategic importance to society. The MILP problems, however, suffer from combinatorial complexity. Because of integer decision variables, as the problem size increases, the number of possible solutions increases super-linearly thereby leading to a drastic increase in the computational effort. To efficiently solve MILP problems, a “price-based” decomposition and coordination approach is developed to exploit 1. the super-linear reduction of complexity upon the decomposition and 2. the geometric convergence potential inherent to Polyak’s stepsizing formula for the fastest coordination possible to obtain near-optimal solutions in a computationally efficient manner. Unlike all previous methods to set stepsizes heuristically by adjusting hyperparameters, the key novel way to obtain stepsizes is purely decision-based: a novel “auxiliary” constraint satisfaction problem is solved, from which the appropriate stepsizes are inferred. Testing results for large-scale Generalized Assignment Problems demonstrate that for the majority of instances, certifiably optimal solutions are obtained. For stochastic job-shop scheduling as well as for pharmaceutical scheduling, computational results demonstrate the two orders of magnitude speedup as compared to Branch-and-Cut. The new method has a major impact on the efficient resolution of complex Mixed-Integer Programming problems arising within a variety of scientific fields.

The associated systems are created by interconnecting I smaller subsystems, each having its own objective and a set of constraints. The subsystem interconnection is modeled through the use of system-wide coupling constraints. Accordingly, the MILP problems are frequently formulated in terms of cost components associated with each subsystem with the corresponding objective functions being additive as such: Furthermore, coupling constraints are additive in terms of I subsystems: The primal problem (1), (2) is assumed to be feasible and the feasible region is assumed to be bounded and finite. The MILP problems modeling the above systems are referred to as separable.
Because of the discrete decisions, however, MILP problems are known to be NP-hard and are prone to the curse of combinatorial complexity. As the size of a problem increases, the associated number of combinations of possible www.nature.com/scientificreports/ solutions (hence the term "combinatorial") increases super-linearly (e.g., exponentially) thereby making problems of practical sizes difficult to solve to optimality; even near-optimal solutions are frequently difficult to obtain. A beacon of hope to resolve combinatorial difficulties lies through the exploitation of separability through the dual "price-based" decomposition and coordination Lagrangian Relaxation technique. After the relaxation of coupling constraints (2), the coordination of subproblems amounts to the maximization of a concave nonsmooth dual function: where Here L(x, y, is the Lagrangian function. The Lagrangian multipliers ("dual" variables) are the decision variables with respect to the dual problem (3), and it is assumed that the set of optimal solutions is not empty. The minimization within (4) with respect to {x, y} is referred to as the "relaxed problem. " While the sizes of the primal and the relaxed problems are the same in terms of the number of discrete variables, the main advantage of Lagrangian Relaxation is the exploitation of the reduction of the combinatorial complexity upon decomposition into subproblems. Accordingly, the number of discrete decision variables within the primal problem is n = I i=1 n x i , so the worst-case complexity of solving the primal problems is O(e I i=1 n x i ) . By the same token, the worst-case complexity required to solve the following subproblem is O(e n x i ) . The decomposition "reverses" the combinatorial complexity thereby exponentially reducing the effort. The decomposition, therefore, offers a viable potential to improve the operations of existing systems as well as to scale up the size of the systems to support their efficient operations.
While decomposition efficiently reduces the combinatorial complexity, the coordination aspect of the method to efficiently obtain the optimal "prices" (Lagrangian multipliers) has been the subject of an intense research debate for decades because of the fundamental difficulties of non-smooth optimization. Namely, because of the presence of integer variables x, the dual function (3) is non-smooth comprised of flat convex polygonal facets (each corresponding to a particular solution to the relaxed problem within (4)) intersecting at linear ridges along which the dual function q( ) is non-differentiable; in particular, q( ) is not differentiable at * thereby ruling out the possibility of using necessary and sufficient conditions for the extremum. As a result of the non-differentiability of q( ) , subgradient multiplier-updating directions, however, are non-ascending directions thereby leading to a decrease of dual values; subgradient directions may also change drastically thereby resulting in zigzagging of Lagrangian multipliers (see Fig. 1 for illustrations) and slow convergence as a result.  www.nature.com/scientificreports/ Traditional methods to maximize q( ) rely upon iterative updates of Lagrangian multipliers by taking a series of steps s k along subgradient g(x k , y k ) directions as: is a an optimal solution to the relaxed problem (4) with multipliers equal to k . Within the Lagrangian Relaxation framework, subgradients are defined as levels of constraint violations if present, can be handled by converting into equality constraints by introducing non-negative real-valued slack variables z The multipliers are subsequently projected onto the positive orthant delineated by restrictions ≥ 0.
Because of the lack of differentiability of q( ) , notably, at the optimum * , the stepsize selection plays an important role to guarantee convergence to the optimum as well as for the success of the overall Lagrangian Relaxation methodology for solving MILP problems.
One of the earlier papers on the optimization of non-smooth convex functions, with q( ) being its member, though irrespective of Lagrangian Relaxation, is Polyak's seminal work 22 . Intending to achieve the geometric (also referred to as "linear") rate of convergence so that � k − * � is monotonically decreasing, Polyak proposed the stepsizing formula, which in terms of the problem under consideration takes the following form: Within (7) and thereafter in the paper the standard Euclidean norm is used.
Subgradient directions, however, 1. are generally difficult to obtain computationally when the number of subproblems (5) to be solved is large, and 2. change drastically thereby resulting in zigzagging of Lagrangian multipliers and slow convergence. Moreover, 3. stepsizes (7) cannot be set due to the lack of the knowledge about the optimal dual value q( * ).
To overcome the first two of the difficulties above, the Surrogate Subgradient method was developed by 23 whereby the exact optimality of the relaxed problem (or even subproblems) is not required. As long as the following "surrogate optimality condition" is satisfied: the multipliers can be updated by using the following version of the Polyak's formula and convergence to * is guaranteed. Here "tilde" is used to distinguish optimal solutions {x k , y k } to the relaxed problem from the solutions {x k ,ỹ k } that satisfy the "surrogate optimality condition" (8). Unlike that in Polyak's formula, parameter γ is less than 1 to guarantee that q( * ) > L(x k ,ỹ k , k ) so that the stepsizing formula (9) is well-defined, as proved by Zhao et al. 23 . Once {x k ,ỹ k } are obtained, multipliers are updated by using the same formula as in (6) with stepsizes from (9) and "surrogate subgradient" multiplier-updating directions g(x k ,ỹ k ) used in place of subgradient directions g(x k , y k ) . Besides reducing the computational effort owing to (8), the concomitant reduction of multiplier zigzagging has also been observed. The main difficulty is the lack of knowledge about q( * ) . As a result, the geometric/linear convergence of the method (or any convergence at all) is highly questionable in practice. Nevertheless, the underlying geometric convergence principle behind the formula (8) is promising and will be exploited in "Results" section.
One of the first attempts to overcome the difficulty associated with the unavailability of the optimal [dual] value is the Subgradient-Level method developed by Goffin and Kiwiel 24 by adaptively adjusting a "level" estimate based on the detection of "sufficient descent" of the [dual] function and "oscillation" of [dual] solutions. In a nutshell, a "level" estimate is set as q k lev = q k j rec + δ j with q k rec being the best dual value ("record objective value") obtained up to an iteration k, and δ j is an adjustable parameter with j denoting the j th update of q k lev . Every time oscillations of multipliers are detected, δ j is reduced by half. In doing so, stepsizes appropriately decrease, q k lev increases (for maximization of non-smooth functions such as (3)) and the process continues until δ j → 0 and q k lev → q( * ).
To improve convergence, rather than updating all the multipliers "at once, " within the Incremental Subgradient methods 25 , multipliers are updated "incrementally. " Convergence results of the Subgradient-Level method 24 have been extended for the Incremental Subgradient methods.
Within the Surrogate Lagrangian Relaxation (SLR) method 26 , the computational effort is reduced along the lines of the Surrogate Subgradient method 23 discussed above, that is, by solving one of a few subproblems at a time. To guarantee convergence, within SLR, distances between multipliers at consecutive iterations are required to decrease through a specially-constructed contraction mapping until convergence. As demonstrated by Bragin et al. 26 , the SLR method converges faster as compared to the above-mentioned Subgradient-Level method 24 and the Incremental Subgradient methods 25,27 for non-smooth optimization. Unlike the Subgradient-Level and Incremental Subgradient methods 25,27 , the SLR method does not require obtaining dual values to set stepsizes, which further reduces the effort. Aiming to simultaneously guarantee convergence while ensuring fast www.nature.com/scientificreports/ reduction of constraint violations and preserving the linearity, the Surrogate Absolute-Value Lagrangian Relaxation (SAVLR) method 28 was developed to penalize constraint violations by using l 1 "absolute-value" penalty terms. The above methods are reviewed in more detail in Supplementry Information Section. Because of the presence of the integer variables, there is the so-called the duality gap, which means that even at convergence, q( * ) is generally less than the optimal cost of the original problem (1), (2). To obtain a feasible solution to (1), (2), the subproblem solutions when put together may not satisfy all the relaxed constraints. Therefore, to solve corresponding MILP problems, heuristics are inevitable and are used to perturb subproblem solutions. The important remark here is that the closer the multipliers are to the optimum, generally, the closer the subproblem solutions are to the global optimum of the original problem, and the easier it is to obtain feasible solutions through heuristics. Therefore, having fast convergence in the dual space to maximize the dual function (3) is of paramount importance for the overall success of the method. Specific heuristics will be discussed at the end of the "Results" section.

Results
Surrogate "Level-Based" Lagrangian Relaxation. In this subsection, a novel Surrogate "Level-Based" Lagrangian Relaxation (SLBLR) method is developed to determine "level" estimates of q( * ) within the Polyak's stepsizing formula (9) for fast convergence of multipliers when optimizing the dual function (3). Since the knowledge of q( * ) is generally unavailable, over-estimates of the optimal dual value, if used in place of q( * ) within the formula (9), may lead to the oscillation of multipliers and to the divergence. Rather than using heuristic "oscillation detection" of multipliers used to adjust "level" values 24 , the key of SLBLR is the decision-based "divergence detection" of multipliers based on a novel auxiliary "multiplier-divergence-detection" constraint satisfaction problem.
"Multiplier-Divergence-Detection" problem to obtain the estimate of q( * ). The premise behind the multiplierdivergence detection is the rendition of the result due Zhao et al. 23 :

Theorem 1 Under the stepsizing formula
such that {x k ,ỹ k } satisfy the multipliers move closer to optimal multipliers * iteration by iteration: The following Corollary and Theorem 2 are the main key results of this paper.

Corollary 1 If then
Theorem 2 If the following auxiliary "multiplier-divergence-detection" feasibility problem (with being a continuous decision variable: ∈ R m ) admits no feasible solution with respect to for some k j and n j , then ∃ κ ∈ [k j , k j + n j ] such that Proof Assume the contrary: ∀κ ∈ [k j , k j + n j ] the following holds: By Theorem 1, multipliers approach * , therefore, the "multiplier-divergence-detection" problem admits at least one feasible solution - * . Contradiction.  (16) it follows that ∃ q κ,j such that q κ,j > q( * ) and the following holds: The equation (18) can equivalently be rewritten as: Therefore, A brief yet important discussion is in order here. The overestimate q j of the dual value q( * ) is the sought-for "level" value after the j th update (the j th time the problem (15) is infeasible). Unlike previous methods, which require heuristic hyperparameter adjustments to set level values, within SLBLR, level values are obtained by using the decision-based principle per (15) precisely when divergence is detected without any guesswork. In a sense, SLBLR is hyperparameter-adjustment-free. Specifically, neither "multiplier-divergence-detection" problem (15), nor the computations within (18)-(20) requires hyperparameter adjustment. Following Nedić and Bertsekas 27 , the parameter γ will be chosen as a fixed value γ = 1 I , which is the inverse of the number of subproblems and will not require further adjustments.
Note that (15) simplifies to an LP constraint satisfaction problem. For example, after squaring both sides of the first inequality � − k j +1 � ≤ � − k j � within (15), after using the binomial expansion, and canceling � − k j � 2 from both sides, the inequality simplifies to To speed up convergence, a hyperparameter ζ < 1 is introduced to reduce stepsizes as follows: Subsequently after iteration k j+1 , the problem (15) is sequentially solved again by adding one inequality per multiplier-updating iteration until iteration k j+1 + n j+1 − 1 is reached for some n j+1 so that (15) is infeasible. Then, stepsize is updated by using q j+1 per (21) and is used to update multipliers until the next time it is updated to q j+2 when the "multiplier-divergence-detection" problem is infeasible again, and the process repeats. Per (21), SLBLR requires hyperparameter ζ , yet, it is set before the algorithm is run and subsequently is not adjusted (see "Numerical testing" section for empirical demonstration of the robustness of the method with respect to the choice of hyperparameter ζ).
To summarize the advantage of SLBLR, hyperparameter adjustment is not needed. The guesswork of when to adjust the level-value, and by how much is obviated -after (15) is infeasible, the level value is formulaically recalculated.
On improvement of convergence. To speed up the acceleration of the multiplier-divergence detection through the "multiplier-divergence-detection" problem, (15) is modified, albeit heuristically, in the following way: Unlike the problem (15), the problem (22) no longer simplifies to an LP problem. Nevertheless, the system of inequalities delineate the convex region and can still be handled by commercial software.
Discussion of (22). Equation (22) is developed based on the following principles: 1. Rather than detecting divergence per (15), convergence with a rate slower than √ 1 − 2 · ν · s is detected. This will lead to a faster adjustment of the level values. While the level value may no longer be guaranteed to be the upper bound to q( * ) , the merit of the above scheme will be empirically justified in the "Numerical testing" section. 2. While the rate of convergence is unknown, in the "worst-case" scenario √ 1 − 2 · ν · s is upper bounded by 1 with ν = 0 , thereby reducing (22) to (15). The estimation of √ 1 − 2 · ν · s is thus much easier than the previously used estimations of q( * ) (as in Subgradient-Level and Incremental Subgradient approaches). 3. As the stepsize approaches zero, √ 1 − 2 · ν · s approaches the value of 1 regardless of the value of ν , once again reducing (22) to (15).
update multipliers per (6) by using g(x k ,ỹ k ) as end if 9: i ← i + 1 10: < ε then search for feasible solutions x f eas , y feas that satisfy (2) to obtain a feasible cost There are three things to note here. 1. Steps in lines 15-16 are optional since other criteria can be used such as the number of iterations or the CPU time; 2. The value of q( k ) is still needed (line 1) to obtain a valid lower bound. To obtain q( k ) , all subproblems are solved optimally for a given value of multipliers k . The frequency of the search for the value q( k ) is determined based on criteria as stated in point 1 above; 3. The search for feasible solutions is explained below.
Search for feasible solutions. Due to non-convexities caused by discrete variables, the relaxed constraints are generally not satisfied through coordination, even at convergence. Heuristics are thus inevitable, yet, they are the last step of the feasible-solution search procedure. Throughout all examples considered, following 28 (as discussed in Supplementary Information), l 1 -absolute-value penalties penalizing constraint violations are considered. After the total constraint violation reaches a small threshold value, a few subproblem solutions obtained by the Lagrangian Relaxation method are perturbed, e.g., see heuristics within accompanying CPLEX codes within 28 to automatically select which subproblem solutions are to be adjusted to eliminate the constraint violation to obtain a solution feasible with respect to the overall problem. Numerical Testing. In this subsection, a series of examples are considered to illustrate different aspects of the SLBLR method. In "Demonstration of convergence of multipliers based on a small example with known optimal multipliers" section, a small example with known corresponding optimal Lagrangian multipliers is considered to test the new method as well as to compare how fast Lagrangian multipliers approach their optimal values as compared to Surrogate Lagrangian Relaxation 26 and to Incremental Subgradient 25 methods. In "Generalized Assignment Problems" section, large-scale instances of generalized assignment problems (GAPs) of types D and E with 20, 40, and 80 machines and 1600 jobs from the OR-library (https:// www-or. amp.i. kyoto-u. ac. jp/ membe rs/ yagiu ra/ gap/) are considered to demonstrate efficiency, scalability, robustness, and competitiveness of the method with respect to the best results available thus far in the literature. In "Stochastic job-shop scheduling with the considerationof scrap and rework" section, a stochastic version of a job-shop scheduling problem instance with 127 jobs and 19 machines based on Hoitomt et al. 29 is tested. In "Multi-stage pharmaceutical scheduling" section, two instances of pharmaceutical scheduling with 30 and 60 product orders, 17 processing units, and 6 stages based on Kopanos et al. 13 are tested.
For "Demonstration of convergence of multipliers based on a small example with known optimal multipliers" section and "Generalized Assignment Problems" section, SLBLR is implemented within CPLEX 12.10 by using a Dell Precision laptop Intel(R) Xeon(R) E-2286M CPU @ 2.40GHz with 16 cores and installed memory (RAM) of 32.0 GB. For "Stochastic job-shop scheduling with the considerationof scrap and rework" section and "Multi-stage pharmaceutical scheduling" section, SLBLR is implemented within CPLEX 12.10 by using a server Intel(R) Xeon(R) Gold 6248R CPU @ 3.00GHz with 48 cores and installed memory (RAM) of 192.0 GB.

Demonstration of convergence of multipliers based on a small example with known optimal multipliers.
To demonstrate convergence of multipliers, consider the following example (due Bragin et al. 30  www.nature.com/scientificreports/ As proved by Bragin et al. 30 , the optimal dual solutions are * 1 = 0.6 and * 2 = 0. Inequality constraints are converted to equality constraints after introducing slack variables. In Fig. 2, the decrease of the corresponding distances from current multipliers to the optimal multipliers ( � k − * � ) is shown, and the SLBLR method is compared with the Incremental Subgradient method 25 and the Surrogate Lagrangian Relaxation method 26 .
Within the SLBLR method, the equation (15) is used to detect divergence, and ζ = 1 2 is used to set stepsizes within (21). In essence, only one hyperparameter was required, which has a quite simple explanation -"when the stepsize is 'too large, ' cut the stepsize in half. " As demonstrated in Fig. 2, the SLBLR method converges fast with � k − * � decreasing roughly along a straight line on a log-scale graph suggesting that the rate of convergence is likely linear as expected.
As for the Incremental Subgradient method, two hyperparameters are required: R and δ (corresponding values used are shown in parentheses in the legend of Fig. 2 (left)). A trial-and-error analysis indicated that "acceptable" values are R = 0.25 and δ = 24. Increasing or decreasing R to 0.5 and 0.125, respectively, do not lead to improvements. Likewise, increasing or decreasing δ to 48 and 12, respectively, do not lead to improvements as well. "Plateau" regions in the figure are caused by the fact that as stepsizes get smaller, a larger number of iterations is required for multipliers to travel the predetermined distance R; during these iterations, stepsizes are not updated and multipliers may oscillate around a neighborhood of the optimum without getting closer. While the above difficulty can be alleviated and convergence can be improved by hyperparameters τ , β , and R l as reviewed in Supplementary Information, however, an even larger number of hyperparameters would be required.
As for the Surrogate Lagrangian Relaxation method, several pairs of hyperparameters (M and r) have been used as well (corresponding values used are shown in parentheses in the legend of Fig. 2 (right)), yet, the performance of Surrogate Lagrangian Relaxaton does not exceed the performance of the SLBLR method.
Herein lies the advantage of the novel SLBLR method: the decision-based principle behind computing the "level" values. This is in contrast to the problem-dependent choice of hyperparameters R and δ within the Subgradient-Level 24 and Incremental Subgradient 25 methods, and the choice of M and r within Surrogate Lagrangian Relaxation 26,28 (see "Introduction" section and Supplementry Information for more detail).
Even after obtaining "appropriate" values of the aforementioned hyperparameters by using a trial-anderror procedure that entails effort, results obtained by Surrogate Lagrangian Relaxation 26 and the Incremental Subgradient method 25 do not match or beat those obtained by the SLBLR method. The specific reasons are 1. Heuristic adjustments of the "level" values are required 24,25 based on multiplier "oscillation detection" or "significant descent" (for minimization of non-smooth functions). However, these rules do not detect whether multipliers "start diverging. " Moreover, oscillation of multipliers is a natural phenomenon when optimizing nonsmooth functions as discussed in "Introduction" section since multipliers may zigzag/oscillate across ridges of the function, so the multiplier "oscillation detection" may not necessarily warrant the adjustment of level values. On the other hand, multiplier "oscillation" is detected by checking whether multipliers traveled a (heuristically) predetermined distance R, hence, the divergence of multipliers can go undetected for a significant number of iterations (hence, the "plateau" regions shown in Fig. 2 (left)), depending on the value of R. To the best of the (23) min www.nature.com/scientificreports/ authors' knowledge, the subgradient-and surrogate-subgradient-based methods using Polyak's stepsizes with the intention of achieving the geometric/linear convergence rate either require q( * ) , which is unavailable, or require multipliers to travel infinite distance to guarantee convergence to the optimum * 24 . 2. While SLR avoids the need to estimate q( * ) , the geometric/linear convergence is only possible outside of a neighborhood of * 26 . Precisely for this reason, the convergence of multipliers within SLR with the corresponding stepsizing parameters M = 30 and r = 0.01 (as shown in Fig. 2 (right)) appears to follow closely convergence within SLBLR up until iteration 50, after which the improvement tapers off.
Generalized assignment problems. To demonstrate the computational capability of the new method as well as to determine appropriate values for key hyperparameters ζ and ν while using standard benchmark instances, large-scale instances of GAPs are considered (formulation is available in subsection 4.2 of Supplementary Information). We consider 20, 40, and 80 machines with 1600 jobs (https://www-or.amp.i.kyoto-u.ac.jp/members/ yagiura/gap/).
To determine values for ζ within (21) and ν within (22) to be used throughout the examples, several values are tested using GAP instance d201600. In Table 1, with fixed values of ν = 2 and s 0 = 0.02 , the best result (both in terms of the cost and the CPU time) is obtained with ζ = 1/1.5 . With the value of ζ = 1/4, the stepsize decreases "too fast" thereby leading to a larger number of iterations and a much-increased CPU time as a result. Likewise, in Table 2 with fixed values of ζ = 1/1.5 and s 0 = 0.02 , it is demonstrated that the best result (both in terms of the cost and the CPU time) is obtained with ν = 2 . Empirical evidence here suggests that the method is stable for other values of ν. The robustness with respect to initial stepsizes ( s 0 ) is tested and the results are demonstrated in Table 3. Multipliers are initialized by using LP dual solutions. The method's performance is appreciably stable for the given range of initial stepsizes used (Table 3). SLBLR is robust with respect to initial multipliers 0 (Table 4). For this purpose, the multipliers are initialized randomly by using the uniform distribution U [90,110]. For the testing, the initial stepsize s 0 = 0.02 was used. As evidenced from Table 4, the method's performance is stable, exhibiting only a slight degradation of solution accuracy and an increase of the CPU time as compared to the case with multipliers initialized by using LP dual solutions.
To test the robustness as well as scalability of the method across several large-scale GAP instances, six instances d201600, d401600, d801600, e201600, e401600, and e801600 are considered. SLBLR is compared with Depth-First Lagrangian Branch-and-Bound method (DFLBnB) 31 , Column Generation 32 , and Very Large Scale Neighborhood Search (VLSNS) 33 , which to the best of the authors' knowledge are the best methods for at least one of the above instances. For completeness, a comparison against Surrogate Absolute-Value Lagrangian Relaxation (SAVLR) 28 , which is an improved version of Surrogate Lagrangian Relaxation (SLR) 26 , is also performed. The latter SLR method 26 has been previously demonstrated to be advantageous against other non-smooth optimization  www.nature.com/scientificreports/ methods as explained in "Introduction" section. Table 5 presents feasible costs and times (in seconds) for each method. The advantage of SLBLR is the ability to obtain optimal results across a wider range of GAP instances as compared to other methods. Even though the comparison in terms of the CPU time is not entirely fair, feasiblecost-wise, SLBLR decisively beats previous methods. For the d201600 instance, the results obtained by SLBLR and the Column Generation method 32 are comparable. For instance d401600, SLBLR obtains a better feasible solution and for instance d801600, the advantage over the existing methods is even more pronounced.
To the best of the authors' knowledge, no other reported method obtained optimal results for instances d401600 and d801600. SLBLR outperforms SAVLR 28 as well, thereby demonstrating that the fast convergence offered by the novel "level-based" stepsizing, with other things being equal, translates into better results as compared to those obtained by SAVLR, which employs the "contraction mapping" stepsizing 28 . Lastly, the methods developed in 31-33 specifically target GAPs, whereas the SLBLR method developed in this paper has broader applicability.
Stochastic job-shop scheduling with the consideration of scrap and rework. To demonstrate the computational capability of the method to solve large-scale stochastic MILP problems, a job-shop scheduling problem is considered. Within a job shop, each job requires a specific sequence of operations and the processing time for each operation. Operations are performed by a set of eligible machines. To avoid late shipments, expected tardiness is minimized. Limited machine capacity brings a layer of difficulty since multiple "individual-job" subproblems are considered together competing for limited resources (machines). Another difficulty arises because of uncertainties, including processing times [34][35][36][37][38][39] and scrap [40][41][42] . Re-manufacturing of one part may affect and disrupt the overall schedule within the entire job shop, thereby leading to unexpectedly high delays in production.
In this paper, we modified data from the paper by Hoitomt et al. 29 by modifying several jobs by increasing the number of operations (e.g., from 1 to 6) and decreasing the capacities of a few machines; the data are in Tables S1 and S2. The stochastic version of the problem with the consideration of scrap and rework is available within the manuscript by Bragin et al. 42 . With these changes, the running time of CPLEX spans multiple days as demonstrated in Fig. 3. In contrast, within the new method, a solution of the same quality as that obtained by CPLEX, is obtained within roughly 1 hour of CPU time. The new method is operationalized by relaxing machine capacity constraints 42 and coordinating resulting job subproblems; at convergence, the beginning times of several jobs are adjusted by a few time periods to remove remaining machine capacity constraint violations.
Multi-stage pharmaceutical scheduling. To demonstrate the capability of the method to solve scheduling problems complicated by the presence of sequence-dependent setup times, a multi-stage pharmaceutical scheduling problem proposed by Kopanos et al. 13 is considered. Setup times vary based on the sequencing of products on each unit (machine). Scheduling in this context is combinatorial in the number of product orders, units, and stages. The new method is operationalized by relaxing constraints that couple individual processing units, Table 4. Robustness results for instance d201600 with respect to initial multipliers 0 . The best feasible cost values obtained are in bold.  Table 5. Comparison against the best results currently available. * The optimality is certified by the LP optimal values, which are 97105 and 97034 for instances d401600 and d801600, respectively. * * The optimality is certified through the lower bound results of, i.e., Posta et al. 31 . − † Not solved to optimality within 24 hours and not reported within the original paper by Posta et al. 31 . − These instances were not considered within the papers by Sadykov et al. 32 and Bragin et al. 28 Fig. 4.

Case number Feasible cost
With a relatively small number of product orders, 30, an optimal solution with a feasible cost of 54.97 was found by CPLEX within 1057.78 seconds. The optimality is verified by running CPLEX until the gap is 0%; it took 171993.27 seconds to verify the optimality. SLBLR takes a slightly longer time to obtain the same solution -1647.35 seconds (Fig. 4 (left)). In contrast, with 60 product orders, CPLEX no longer obtains good solutions in a computationally efficient manner; a solution with a feasible cost of 55.98 is obtained after 1,000,000 seconds. Within SLBLR, a solution with a feasible cost of 55.69 is obtained within 1978.04 seconds. This constitutes more than two orders of magnitude of improvement over CPLEX as demonstrated in Fig. 4 (right; log scale). When  www.nature.com/scientificreports/ doubling the number of products, CPLEX's performance is drastically deteriorated, while the performance of SLBLR is scalable.

Discussion
This paper develops a novel MILP solution methodology based on the Lagrangian Relaxation method. Salient features of the novel SLBLR method, inherited from the previous versions of Lagrangian Relaxation, are: 1. reduction of the computational effort required to obtain Lagrangian-multiplier-updating directions and 2. alleviation of zigzagging of multipliers. The key novel feature of the method, which the authors believe gives SLBLR the decisive advantage, is the innovative exploitation of the underlying geometric-convergence potential inherent to Polyak's stepsizing formula without the heuristic adjustment of hyperparameters for the estimate of q( * ) -the associated "level" values are determined purely through the simple auxiliary "multiplier-divergence-detection" constraint satisfaction problem. Through testing, it is discovered that SLBLR is robust with respect to the choice of initial stepsizes and multipliers, computationally efficient, competitive, and general. Several problems from diverse disciplines are tested and the superiority of SLBLR is demonstrated. While "separable" MILP problems are considered, no particular problem characteristics such as linearity or separability have been used to obtain "level" values, and thus SLBLR has the potential to solve a broad class of MIP problems.