EGD-ILS: Hybrid Energy Gradient Descent Optimization for Job Shop Scheduling

1. Introduction

The JSSP is a central problem in operations research and production management. It consists of assigning a set of n jobs on a set of m machines. Each job $J_j$ consists of m operations $O_{jo}$ that must be processed in a given order. An operation $O_{jo}$ has a processing time $p_{jo}$ and it has a predefined machine $M_{jo}$. Our goal is to achieve the minimum makespan. Precedence and resource constraints must be satisfied.

This problem is NP-hard. Exact approaches such as Integer Linear Programming or Constraint Programming fail when the size exceeds 10x10. Metaheuristics such as Genetic Algorithms, Simulated Annealing, or Particle Swarm Optimization provide good results but require fine parameter tuning and spend much time in infeasible regions of the search space [1].

Since the late 1980s, Hopfield networks were applied to JSSP by encoding each decision as a binary variable $V_{jot}$ indicating if job j is processed on operation o at time t [2]. This leads to $O\left(nmT\right)$ neurons with T is the number of time steps. With Deep Learning, GNN (Graph Neural Networks) and RL (Reinforcement Learning) methods achieved strong results but require GPU and large datasets [3]. Recently, hybrid approaches have emerged to combine the advantages of different methods. The approach involves breaking down the problem into two parts: 1. Find a feasible solution, 2. Optimize the objective function [4].

We propose EGD-ILS, a hybrid metaheuristics that uses the energy-based philosophy but optimize directly the continuous vector of start times $S\in \mathbb{R}{+}^{n\times m}$. A differentiable energy $E\left(S\right)$ measures the magnitude of violations. Minimizing $E\left(S\right)$ via gradient descent repairs an infeasible solution until an admissible schedule is obtained. Iterated Local Search explores admissible solutions to find optimal or near optimal solution. Our main contributions are threefold:

1.

We prove mathematically that by fixing the operation sequence, the feasibility search becomes a convex optimization problem that converges to E=0.

2.

We propose a simple and efficient architecture that combines the speed of Gradient Descent for feasibility and the power of Iterated Local Search for optimization.

3.

Experimental validation on FT and LA benchmarks.

4.

Our work is the first work that provides a theoretical convergence guarantee for the feasibility phase of a hybrid JSSP solver based on gradient methods.

The rest of the paper is laid out this way. Section 2 presents the literature review. Section 3 presents the EGD-ILS methodology with the convergence proof. Section 4 discus experimental results. Section 5 concludes.

2. Review and Related Work

The use of neural-inspired energy functions for scheduling dates back to the late 1980. In a pioneering study, Foo and Takefuji [2] formulated JSSP as a minimization problem on a Hopfield network. Their objective combined penalty terms for constraint violations with a term related to the makespan. A key drawback reported was the difficulty in balancing the different penalty weights, which often led to poor convergence behaviour [4]. Van den Bout [5] introduced Mean Field Annealing but this approach suffers from the problem of combinatorial explosion and local minima..

Recent research has focused on energy formulations using modern optimization strategies. for example is the method developed by Peng and others, using an augmented Lagrangian framework which incorporates any constraint violations within the Lagrangian. This approach employed a gradient technique to minimize Lagrangian and then applied a heuristic to achieve feasibility. This approach works well for medium instances only, but it still depends on choice of the correct penalty [6].

To enhance the robustness level, Wang and Zheng proposed a metaheuristic in two phases. In the first one, the simulated annealing is utilized to minimize the compound function of the makespan and the violations of constraints. However, once again, some admissible solution is found, the second phase is based on Tabu search and aimed to minimize the makespan [3,7,8].

In contrast to the methods above, the proposed EGD-ILS does not rely on stochastic search or data-driven training in the first phase. Instead, we show that by committing to a specific operation order upfront, the feasibility problem becomes convex. This allows us to use deterministic Gradient Descent with a proven convergence property. To our knowledge, this explicit convexity argument and its use for guaranteed feasibility in JSSP have not been explored before.

3. Methodology: The EGD-ILS Approach

3.1. Problem Formulation

We developed a mathematical statement of the JSSP problem [9]. Without loss of generality, Let’s use the instance Table 1 as example highlight the mathematical model used.

$S_{ij}$ is the starting time of operation $O_{ij}$. $M_{ij}$ is the resource used. $P_{ij}$ is the processing time . H : Sum of all processing times $H={\sum }_{i=1}^n{\sum }_{j=1}^m{P}_{ij}$.

We introduce the Boolean variable:

$Y_{ip}=1if{S}_{ij}\le S_{pw}and{M}_{ij}=M_{pw}$       $Y_{ip}=0if{S}_{pw}\le S_{ij}and{M}_{ij}=M_{pw}$

The goal is: minimize $C_{max}=max\{S_{ij}+P_{ij}\}$    subject to :

▪ Starting Time constraints : $n\times m$ equations

$$S_{11}\ge 0;S_{12}-8\ge 0;S_{13}-20\ge 0S_{21}\ge 0;S_{22}-11\ge 0;S_{23}-16\ge 0$$

(1)

▪ Precedence Constraints : $n\times (m-1)$ equations

$$S_{12}-S_{11}-8\ge 0;S_{13}-S_{12}-12\ge 0S_{22}-S_{21}-11\ge 0;S_{23}-S_{22}-5\ge 0$$

(2)

▪ Machine Constraints : $n\times m\times (n-1)$ equations

$$S_{22}-S_{11}+H(1-Y_{12})-8\ge 0;S_{11}-S_{22}+H\times Y_{12}-5\ge 0$$

(3)

$$S_{23}-S_{12}+H(1-Y_{12})-12\ge 0;S_{12}-S_{23}+H\times Y_{12}-14\ge 0$$

(4)

$$S_{21}-S_{13}+H(1-Y_{12})-3\ge 0;S_{13}-S_{21}+H\times Y_{12}-11\ge 0$$

(5)

3.2. Phase 1: Feasible Construction by Gradient Descent with Fixed Sequence

3.2.1. Starting Energy

$$\begin{array}{ccc}\hfill S=\left(\begin{array}{c}{S}_{11}\\ S_{12}\\ S_{13}\\ S_{21}\\ S_{22}\\ S_{23}\end{array}\right)& W_s=I_6=\left(\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\hfill & \hfill B_s=\left(\begin{array}{c}0\\ -8\\ -20\\ 0\\ -11\\ -16\end{array}\right)\end{array}$$

(6)

$$Z_s=W_s\times S+B_s$$

(7)

$Z_s$ must be non negative.

$$A_s=max(0;-W_s\times S-B_s)$$

(8)

We define the starting energy as:

$$E_s\left(S\right)={\displaystyle \sum _{i=1}^{n\times m}}{\left(A_s(i,1)\right)}^2$$

(9)

3.2.2. Precedence Energy

$$\begin{array}{ccc}\hfill S=\left(\begin{array}{c}{S}_{11}\\ S_{12}\\ S_{13}\\ S_{21}\\ S_{22}\\ S_{23}\end{array}\right)& W_p=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ -1& 1& 0& 0& 0& 0\\ 0& 0& -1& 1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& -1& 1& 0\\ 0& 0& 0& 0& -1& 1\end{array}\right)\hfill & \hfill B_p=\left(\begin{array}{c}0\\ -8\\ -12\\ 0\\ -11\\ -5\end{array}\right)\end{array}$$

(10)

$$Z_p=W_p\times S+B_p$$

(11)

$Z_p$ must be non negative.

$$A_p=max(0;-W_p\times S-B_p)$$

(12)

We define the precedence energy as:

$$E_p\left(S\right)={\displaystyle \sum _{i=1}^{n\times m}}{\left(A_p(i,1)\right)}^2$$

(13)

3.2.3. Machine Conflict Energy

The key idea of EGD-ILS is to fix the operation sequence in phase 1. This transforms the disjunctive machine constraint into a linear constraint. Without any loss of generality, consider that $\pi =\left\{1;2;2;1;2;1\right\}$ is a fixed operation sequence. Machine Constraints becomes:

$$S_{22}-S_{11}-8\ge 0;S_{11}-S_{22}+H-5\ge 0$$

(14)

$$S_{23}-S_{12}-12\ge 0;S_{12}-S_{23}+H-14\ge 0$$

(15)

$$S_{21}-S_{13}+H-3\ge 0;S_{13}-S_{21}-11\ge 0$$

(16)

$$\begin{array}{ccc}\hfill S=\left(\begin{array}{c}{S}_{11}\\ S_{12}\\ S_{13}\\ S_{21}\\ S_{22}\\ S_{23}\end{array}\right)& W_m=\left(\begin{array}{cccccc}-1& 0& 0& 0& 1& 0\\ 1& 0& 0& 0& -1& 0\\ 0& -1& 0& 0& 0& 1\\ 0& 1& 0& 0& 0& -1\\ 0& 0& -1& 1& 0& 0\\ 0& 0& 1& -1& 0& 0\end{array}\right)\hfill & \hfill B_m=\left(\begin{array}{c}-8\\ H-5\\ -12\\ H-14\\ H-3\\ -11\end{array}\right)\end{array}$$

(17)

$$Z_m=W_m\times S+B_m$$

(18)

$Z_m$ must be non negative.

$$A_m=max(0;-W_m\times S-B_m)$$

(19)

We define the Machine Conflict Energy as:

$$E_m\left(S\right)={\displaystyle \sum _{i=1}^{n\times m}}{\left(A_m(i,1)\right)}^2$$

(20)

3.2.4. Resolution Algorithm

We define the energy function to minimize:

$$E(S,\pi )=E_s\left(S\right)+E_p\left(S\right)+E_m(S,\pi )$$

(21)

Algorithm 1 Energy Descent for JSSP

Input: Instance , $\pi$ , Maxiter , $\alpha$: Learning rate Calculate $W_s,B_s,W_p,B_p,W_m,B_m$ Initialize $S=-B_s$ for $k=1$ to $Maxiter$ do    Compute E(S,$\pi$)    if $E(S,\pi )=0$ then     Return S    else     $[E,g]\leftarrow$ Compute Gradient(S)     $S\leftarrow S-\alpha ·g$     $S\leftarrow max(0,S)$    end if end for Return S

3.2.5. Mathematical Proof of Convergence

Lemma 1.

The energy function $E(S,\pi )$, which was defined above for a fixed order π, is convex and L-smooth.

Proof. 

$E(S,\pi )$ has the form ${\left[max(0,-W\times S-B)\right]}^2$. $g\left(z\right)=max{\left[(0,Z)\right]}^2$ is convex and 2-smooth. $g^{\prime }\left(z\right)=2\left[max(0,Z)\right]$. $|g^{\prime }\left(Z_1\right)-g^{\prime }\left(Z_2\right)|\le 2|Z_1-Z_2|$ . Z is an affine function of S. The composition is convex and L-smooth.

The function $E(S,\pi )$ is a sum of convex and smooth functions. Therefore $E(S,\pi )$ is convex and L-smooth [10].    □

Lemma 2.

For any fixed order π, there exists $S^{*}$ such that $E(S^{*},\pi )=0$.

Proof. 

By definition of the JSSP, for any order $\pi$ there exists at least one feasible schedule. For this schedule, all violations are $\le 0$. As $E\ge 0$, $E=0$ is the global minimum.    □

Theorem 1.

For any step $0<\alpha <2/L$, $S\left(t\right)$ converges to $S^{*}$ such that $E{(}^{*})=0$.

Proof. 

$S(t+1)=max(0,S(t\left)\right)-\alpha \nabla E\left(S\right(t\left)\right)$ From Lemma 1, we can see that E is convex and L-smooth. From Lemma 2, we know that the set of minima is non-empty. Moreover $\nabla E$ is L-Lipschitz. We apply the classical convergence theorem of gradient descent for convex functions. ([11], Sec. 9.1)    □

3.3. Phase 2: Optimization by Iterated Local Search

The ILS explores the set of admissible schedules by using neighbours structure and deciding to transit or not to the candidate solution.

Algorithm 2 Optimization by ILS

Input: Instance ; $\beta$ ; $T_0$ Initialize ${\pi }_0$ ; $\left[S\right]⟵Algorithm1(Instance;{\pi }_0)$ ${\pi }^{\left(best\right)}⟵{\pi }_0$ ;    $\pi ⟵{\pi }_0$ ;    $T⟵T_0$ Compute $C_{max}\left(\pi \right)$ for $k=1$ to K do    Choose $\pi ″$ from neighbours of $N_5\left(\pi \right)$    $\left[S\right]⟵Algorithm1(Instance;\pi ″)$    Compute $C_{max}(\pi ″)$    if $C_{max}(\pi ″)<C_{max}\left(\pi \right)$ then     $\pi$ := $\pi ″$     if $C_{max}\left(\pi \right)<C_{max}\left({\pi }^{\left(best\right)}\right)$ then       ${\pi }^{\left(best\right)}$ := $\pi$     end if    else     $\mathsf{\Delta }=C_{max}(\pi ″)-C_{max}\left(\pi \right)$     if $random\le exp(-\mathsf{\Delta }/T)$ then       $\pi$ := $\pi ″$     end if    end if    T := $\beta \times T$ end for

Denote $N_1$ permutes two operations.$N_2$ is right shift insertion. $N_3$ uses left shift insertion. $N_4$ is inversion blocks between two positions [12]. $N_5$ uses random combination of the four neighbours $N_k\left(\pi \right),k=1,...,4$.

An initial temperature of $T_0$ is given. The parameter $\beta \in [0.5;0.95]$ the speed of convergence.

4. Experimental Results

The RPD is: $RPD={\displaystyle \frac{C_{max}-C_{max}^{*}}{C_{max}^{*}}}\times 100$ where $C_{max}$ is the makespan found by the EGD-ELS; and $C_{max}^{*}$ is the optimal makespan.

The Table 2 and Table 3 present the results on FT and LA benchmarks.

EGD-ILS delivers high-quality solutions on FT and LA benchmarks with reduced computation time. It achieves RPD < 3.9%, solving FT06 and LA05 to LA15 optimally.

5. Conclusions

We proposed a straightforward energy-based approach for JSSP that’s simple and fast. EGD-ILS shows strong performance on FT and LA benchmarks with RPD < 3.9% , solving FT06 and LA05-LA15 optimally with reduced computation time.