Multi-objective Optimization of Energy-Efficient Multi-stage, Multi-level Assembly Job Shop Scheduling

1. Introduction

Since the manufacturing industry is one of the most energy-intensive sectors, how to effectively improve energy efficiency while maintaining production efficiency is currently a major focus within the industry. Energy-efficient scheduling (EES) utilizes limited resources to rationally arrange production and optimize the production process, which is an effective way to effectively reduce the energy consumption in manufacturing systems. Recently, EES has been successfully applied to various scheduling problems in discrete manufacturing industries, yielding positive results.

As an extension of the classic job shop scheduling problem (JSP), the multi-stage, multi-level assembly job shop scheduling problem (MsMlAJSP), which not only involves multiple stages of production from processing to assembly, but also considers complex assemblies composed of multiple jobs with a tree-like multi-level structure consisting of components and subassemblies. It plays an important role in quickly responding to the manufacturing of customized products with multi-variety and small batch. Wang et al. [1] addressed the distributed MsMlAJSP considering maintenance, aiming to minimize the maximum tardiness. Cheng et al. [2] proposed batch transfer and mixed-model assembly for MsMlAJSP with the optimization objectives of minimizing manufacturing and transportation costs. Building on this, Cheng et al. [3] further considered lot streaming with the objective of minimizing total completion time and production cost. All of the above studies have already conducted relevant work on MsMlAJSP, primarily focusing on improving production efficiency and economic benefits. However, as the number of product structural levels and manufacturing stages increases, constraints become more complex, leading to idle time between and within stages, which results in significant idle energy consumption. Therefore, further research is needed on the energy-efficient multi-stage, multi-level assembly job shop scheduling problem (EEMsMlAJSP) for complex multi-level products.

The EEMsMlAJSP is highly complex with a vast solution space, presenting challenges for efficient resolution. Currently, solving the energy-efficient scheduling problem in complex manufacturing systems mainly focuses on energy-efficient strategies and algorithm design. Among these, the block-moving energy-efficient strategy, which makes operations within stages more compact by moving process operations, can effectively reduce idle energy consumption. This is currently an effective method for handling multi-operation and multi-stage complex scheduling problems [4,5,6,7]. There are two main approaches: one is to apply the block-moving strategy only to the solutions in the non-dominated set to further reduce energy consumption [8]; however, this operation may cause the loss of some potential solutions, resulting in a local rather than a global optimum. The other approach is to embed the block-moving strategy into the individual decoding process to avoid losing the global optimal solution [9]; however, such design can lead to a waste of computational resources and severely slow down the search efficiency of algorithm. Therefore, it is an urgent need to design an efficient block-moving strategy that improves solution quality and search efficiency without losing the global optimum.

In terms of algorithm design, the imperialist competitive algorithm (ICA) is an effective evolutionary algorithm framework. It facilitates population evolution through mechanisms such as the assimilation, colonial revolutions, and imperial competition, demonstrating efficient global and local search capabilities. For the past few years, ICA has been successfully applied to solving various multi-objective energy-efficient scheduling problems, such as multi-objective energy-efficient flexible job shop scheduling with energy consumption thresholds [10], energy-efficient hybrid flow shop scheduling [11], energy-efficient flexible job shop scheduling considering transportation and setup times [12], and energy-efficient flexible job shop batch scheduling [13]. For the assimilation in the global search of the algorithm, existing studies mostly adopt the Uniform Crossover (UX) operator, with the probability parameters of uniform crossover preset through experiments. However, applying traditional static parameters to ICA may not ensure optimal search efficiency and solution quality due to the complexity of the EEMsMlAJSP problem. Therefore, it is necessary to design a reasonable dynamic parameter adaptive adjustment mechanism to balance the exploration and exploitation of the algorithm. Reinforcement Learning (RL), such as Q-learning, offers strong decision-making and optimization capabilities. Recent studies have shown that Q-learning can effectively adjust algorithm parameters, enhancing the balance between exploration and exploitation and improving solution quality [14,15,16].

In summary, we propose an improved ICA based on Q-learning (IICA-QL) to solve the EEMsMlAJSP, aiming to the maximum completion time and total energy consumption. First, a mathematical model of the problem is established, and an energy-efficient strategy decoding trigger mechanism is designed to simultaneously enhance solution quality and search efficiency. Second, the IICA-QL for solving the EEMsMlAJSP is proposed. In IICA-QL, two heuristic rules and a random generation method are used to collaboratively construct an initial population with high quality and diversity. A dynamic adaptive adjustment strategy for assimilation operator parameter is implemented based on Q-learning to enhance convergence speed while maintaining population diversity, thus the exploration and exploitation capabilities of the algorithm are better balanced. A hyper-heuristic variable neighborhood search-guided revolution operation is designed to conduct detailed searches on identified high-quality colonies. The competitive phase is replaced with a joint imperialist invasion operation to achieve cooperative co-evolution and information exchange among multiple empires, resulting in more evenly distributed non-dominated solutions. Finally, the effectiveness of the proposed algorithm in solving the EEMsMlAJSP is verified through experimental results and comparison with other algorithms.

The specific structure of the paper is as follows: The problem definition and Mathematical model are established in Section 2. The encoding and decoding design is presented in Section 3. The introduction of the improved imperialist competitive algorithm based on Q-learning (IICA-QL) is described in Section 4. Experimental results and analysis are provided in Section 5. The conclusion to the paper is given in Section 6.

2. Problem Formulation

Table 1. Notations of mathematical model.

Notations	Descriptions
Indices
$z$	Index of products, $z=1,2,\cdots ,N_P$
$j$	Index of jobs, $j=1,2,\cdots ,N_Q$
$i$	Index of processing operations, $i=1,2,\cdots ,N_O$
$a$	Index of assembly operations, $a=1,2,\cdots ,N_A$
$m$	Index of processing machines, $m=1,2,\cdots ,N_M$
$t$	Index of transportation vehicles, $t=1,2,\cdots ,N_T$
$w$	Index of transportation vehicle load transportation, $w=1,2,\cdots ,N_{tw}$
$k$	Index of assembly machines, $k=1,2,\cdots ,N_K$
Parameters
$M$	The set of processing machines, $M=\{M_1,\cdots ,M_m,\cdots ,M_{N_M}\}$
$V$	The set of transportation vehicles, $V=\{V_1,\cdots ,V_t,\cdots ,V_{N_T}\}$
$K$	The set of assembly machine, $K=\{K_1,\cdots ,K_k,\cdots ,K_{N_K}\}$
$P$	The set of products, $P=\{P_1,\cdots ,P_z,\cdots ,P_{N_P}\}$
$N_p$	Number of products
$N_Q$	Number of jobs
$N_O$	Number of processing operations
$N_A$	Number of assembly operations
$N_M$	Number of processing machines
$N_T$	Number of transportation vehicles
$N_K$	Number of assembly machines
$N_O^m,N_A^k$	Number of operations on processing machine m and assembly machine k
$N_{tw}$	Number of load transportations by vehicle t
$G_{max}$	Maximum load capacity of the vehicle
$J_j,J_z$	The job indexed by j, the set of all jobs belonging to$P_z$
$O_j^i,O_j$	The i-th processing operation of $J_j$, the set of all operations of $J_j$
$A_z^a,A_z$	The a-th assembly operation and the set of assembly operations of $P_z$
$M_j^i$	The set of machines available for processing operation $O_j^i$
$P{P}_m,P{P}_k$	Unit power consumption of processing machine m and assembly machine k in working mode
$P{I}_m,P{I}_k$	Unit power consumption of processing machine m and assembly machine k in standby mode
$P{S}_m$	Unit power consumption of processing machine m in setup mode
$P{T}_t$	Unit power consumption of transportation vehicle t
$Di{s}_{P-A}$	Distance between the processing shop and the assembly shop
$V_0,V_t^w$	Speed of transportation vehicles when unloaded and speed of transportation vehicle t during the w-th load transportation
$P({\pi }_i^m),P({\pi }_a^k)$	Processing time of operation ${\pi }_i^m$ and ${\pi }_a^k$
$SP({\pi }_{i-1}^m,{\pi }_i^m)$	Setup time of operation ${\pi }_{i-1}^m$ and ${\pi }_i^m$
$IP({\pi }_{i-1}^m,{\pi }_i^m)$	Idle time between operation ${\pi }_{i-1}^m$ and ${\pi }_i^m$ excluding setup time
$IP({\pi }_{a-1}^k,{\pi }_a^k)$	Idle time between operation ${\pi }_{a-1}^m$ and ${\pi }_a^m$
$RT({\pi }_1^m)$	Release time of job ${\pi }_1^m$
Decision variables
${\pi }^m=[{\pi }_1^m,{\pi }_2^m,\cdots ,{\pi }_{N_O^m}^m]$	Sequence of operations on processing machine m
${\pi }^k=[{\pi }_1^k,{\pi }_2^k,\cdots ,{\pi }_{N_A^k}^k]$	Sequence of operations on assembly machine k
${\pi }^{st1},{\pi }^{st2},{\pi }^{st3}$	Sequence of operations on processing stage, transportation stage, and assembly stage
$ST({\pi }_a^k),CT({\pi }_a^k)$	Start time and completion time of operation ${\pi }_a^k$
$E{C}_R^{st1},E{C}_S^{st1},E{C}_P^{st1},E{C}_I^{st1}$	Release energy consumption, setup energy consumption, processing energy consumption, and idle energy consumption in the processing stage
$E{C}_A^{st3},E{C}_I^{st3}$	assembly energy consumption and idle energy consumption in the assembly stage
$E{C}^{st1},E{C}^{st2},E{C}^{st3}$	Energy consumption in the processing stage, transportation stage, and assembly stage
$C_{\max }$	Maximum completion time of a schedule
TEC	Total energy consumption

Table 1. Notations of mathematical model.

Notations	Descriptions
Indices
$z$	Index of products, $z=1,2,\cdots ,N_P$
$j$	Index of jobs, $j=1,2,\cdots ,N_Q$
$i$	Index of processing operations, $i=1,2,\cdots ,N_O$
$a$	Index of assembly operations, $a=1,2,\cdots ,N_A$
$m$	Index of processing machines, $m=1,2,\cdots ,N_M$
$t$	Index of transportation vehicles, $t=1,2,\cdots ,N_T$
$w$	Index of transportation vehicle load transportation, $w=1,2,\cdots ,N_{tw}$
$k$	Index of assembly machines, $k=1,2,\cdots ,N_K$
Parameters
$M$	The set of processing machines, $M=\{M_1,\cdots ,M_m,\cdots ,M_{N_M}\}$
$V$	The set of transportation vehicles, $V=\{V_1,\cdots ,V_t,\cdots ,V_{N_T}\}$
$K$	The set of assembly machine, $K=\{K_1,\cdots ,K_k,\cdots ,K_{N_K}\}$
$P$	The set of products, $P=\{P_1,\cdots ,P_z,\cdots ,P_{N_P}\}$
$N_p$	Number of products
$N_Q$	Number of jobs
$N_O$	Number of processing operations
$N_A$	Number of assembly operations
$N_M$	Number of processing machines
$N_T$	Number of transportation vehicles
$N_K$	Number of assembly machines
$N_O^m,N_A^k$	Number of operations on processing machine m and assembly machine k
$N_{tw}$	Number of load transportations by vehicle t
$G_{max}$	Maximum load capacity of the vehicle
$J_j,J_z$	The job indexed by j, the set of all jobs belonging to$P_z$
$O_j^i,O_j$	The i-th processing operation of $J_j$, the set of all operations of $J_j$
$A_z^a,A_z$	The a-th assembly operation and the set of assembly operations of $P_z$
$M_j^i$	The set of machines available for processing operation $O_j^i$
$P{P}_m,P{P}_k$	Unit power consumption of processing machine m and assembly machine k in working mode
$P{I}_m,P{I}_k$	Unit power consumption of processing machine m and assembly machine k in standby mode
$P{S}_m$	Unit power consumption of processing machine m in setup mode
$P{T}_t$	Unit power consumption of transportation vehicle t
$Di{s}_{P-A}$	Distance between the processing shop and the assembly shop
$V_0,V_t^w$	Speed of transportation vehicles when unloaded and speed of transportation vehicle t during the w-th load transportation
$P({\pi }_i^m),P({\pi }_a^k)$	Processing time of operation ${\pi }_i^m$ and ${\pi }_a^k$
$SP({\pi }_{i-1}^m,{\pi }_i^m)$	Setup time of operation ${\pi }_{i-1}^m$ and ${\pi }_i^m$
$IP({\pi }_{i-1}^m,{\pi }_i^m)$	Idle time between operation ${\pi }_{i-1}^m$ and ${\pi }_i^m$ excluding setup time
$IP({\pi }_{a-1}^k,{\pi }_a^k)$	Idle time between operation ${\pi }_{a-1}^m$ and ${\pi }_a^m$
$RT({\pi }_1^m)$	Release time of job ${\pi }_1^m$
Decision variables
${\pi }^m=[{\pi }_1^m,{\pi }_2^m,\cdots ,{\pi }_{N_O^m}^m]$	Sequence of operations on processing machine m
${\pi }^k=[{\pi }_1^k,{\pi }_2^k,\cdots ,{\pi }_{N_A^k}^k]$	Sequence of operations on assembly machine k
${\pi }^{st1},{\pi }^{st2},{\pi }^{st3}$	Sequence of operations on processing stage, transportation stage, and assembly stage
$ST({\pi }_a^k),CT({\pi }_a^k)$	Start time and completion time of operation ${\pi }_a^k$
$E{C}_R^{st1},E{C}_S^{st1},E{C}_P^{st1},E{C}_I^{st1}$	Release energy consumption, setup energy consumption, processing energy consumption, and idle energy consumption in the processing stage
$E{C}_A^{st3},E{C}_I^{st3}$	assembly energy consumption and idle energy consumption in the assembly stage
$E{C}^{st1},E{C}^{st2},E{C}^{st3}$	Energy consumption in the processing stage, transportation stage, and assembly stage
$C_{\max }$	Maximum completion time of a schedule
TEC	Total energy consumption

2.1. Problem Definition

The EEMsMlAJSP addressed in this work consists of multiple stages, including processing, transportation, and assembly. The problem is defined as follows. In the processing stage, transportation stage, and assembly stage, there are$N_M$ flexible processing machines $M=\{M_1,\cdots ,M_m,\cdots ,M_{N_M}\}$, $N_T$ transportation vehicles $V=\{V_1,\cdots ,V_t,\cdots ,V_{N_T}\}$ with the same load capacity $G_{max}$, and $N_K$ assembly machines $K=\{K_1,\cdots ,K_k,\cdots ,K_{N_K}\}$, respectively. There is a product set $P=\{P_1,\cdots ,P_z,\cdots ,P_{N_P}\}$ that needs to be manufactured, where each product $P_z\in P$ is composed of a set of various jobs $J_z$ and should undergo an assembly operation set $A_z$ to be assembled into the final product. Each job $J_j\in J_z$ has a release time and consists of a set of processing operations $O_j$. all operations $O_j^i\in O_j$ must be completed before the $J_j$ can be transported in batches by vehicles with load capacity constraint. Once all jobs $J_j$​ belonging to a product $P_z$ have been transported to the assembly stage, the assembly of $P_z$ can begin. In the processing stage, each operation $O_j^i$ has a set of compatible machines $M_j^i$ from which one machine is selected for processing. Once the operation is completed, the next operation for job $J_j$ processed according to the specified sequence. In the transportation stage, a fleet of vehicles $V$ is available to transport jobs from the processing stage to the assembly stage. Each vehicle $V_t\in V$ departs from the processing shop and must return after completing its transportation tasks. The vehicles are subject to capacity constraints and are allowed to make multiple trips, with each trip carrying a batch of components produced in the processing stage. In the assembly stage, once all components belonging to product $P_z$​ are ready, all operations​ $A_z^a\in A_z$ are sequentially completed according to the assembly process.

The optimization objectives of this problem are to simultaneously minimize $C_{\max }$ and TEC. It is a multi-objective problem and encompasses several sub-problems, including the sequencing of processing operations, assignment of processing machines, sequencing of job transportation, assignment of vehicles, sequencing of product assembly, sequencing of assembly operations and assignment of assembly machines. The computational complexity of the mathematical model is $O=((N_O!)\times {N_M}^{N_O}\times (N_Q!)\times {N_T}^{N_Q}\times (N_P!)\times {N_K}^{N_A})$. So as the problem size increases, the complexity grows exponentially. Moreover, as an extension of the well-known NP-hard problem JSP, the EEMsMlAJSP integrates assembly and batch transportation, making it an NP-hard problem as well, and even more complex.

The schematic diagram of the EEMsMlAJSP model is shown in Figure 1. Some assumptions are listed below:

• In the processing and assembly stages, any machine can only process one operation at a time.
• Once started, processing operations, transportation processes, and assembly operations cannot be interrupted.
• All machines and vehicles are available at time zero.
• Transportation time between machines for consecutive processing (or assembly) operations is not considered.
• Vehicles maintain a constant speed during transportation.
• The initial setup time is included in the release time of the workpiece.

2.2. Mathematical Model

The mathematical model for the three stages with two optimization objectives in EEMsMlAJSP is constructed as follows:

The first objective, the maximum completion time, is the time when all products are fully assembled, calculated as follows:

$$C_{\max }=\max \{CT({\pi }_a^k)\},\forall a,k$$

(1)

The second objective, namely total energy consumption, is composed of three parts: total energy consumption in the processing stage, total energy consumption in the transportation stage, and total energy consumption in the assembly stage. The total energy consumption for each stage is calculated as follows:

(1) Processing stage: The total energy consumption in this stage is composed of the energy consumption of the processing machines, which is determined by four modes of machines: release mode, setup mode, working mode, and idle mode. The energy consumption for each mode of the processing machines is constructed as follows:

Release energy consumption: This refers to the energy consumption of machines before the first arrival time of the job. The release energy consumption in the processing stage is calculated as follows:

$$E{C}_R^{st1}={\displaystyle \sum _{m=1}^{N_M}P{P}_m\times RT({\pi }_1^m)}$$

(2)

Setup energy consumption: This energy consumption generated during the preparation process of the machine before entering the working mode. The setup energy consumption in the processing stage is calculated as follows:

$$E{C}_S^{st1}={\displaystyle \sum _{m=1}^{N_M}{\displaystyle \sum _{i=2}^{N_{Om}}P{S}_m\times SP({\pi }_{i-1}^m,{\pi }_i^m)}}$$

(3)

Working energy consumption: The energy consumption generated when the processing machine is in working mode. The working energy consumption in the processing stage is calculated as follows:

$$E{C}_P^{st1}={\displaystyle \sum _{m=1}^{N_M}{\displaystyle \sum _{i=1}^{N_{Om}}P{P}_m\times P({\pi }_i^m)}}$$

(4)

Idle energy consumption: The energy consumption generated during the time between processing two consecutive operations, excluding setup time. The idle energy consumption in the processing stage is calculated as follows:

$$E{C}_I^{st1}={\displaystyle \sum _{m=1}^{N_M}{\displaystyle \sum _{i=2}^{N_{Om}}P{I}_m\times IP({\pi }_{i-1}^m,{\pi }_i^m)}}$$

(5)

The total energy consumption in the processing stage is calculated as follows:

$$E{C}^{st1}=E{C}_R^{st1}+E{C}_S^{st1}+E{C}_P^{st1}+E{C}_I^{st1}$$

(6)

(2) Transportation stage: The total energy consumption in this stage consists of the energy consumption of the transportation vehicles, which is composed of two parts: loaded energy consumption and unloaded energy consumption. The total energy consumption in the transportation stage is calculated as follows:

$$E{C}^{st2}={\displaystyle \sum _{t=1}^{N_T}{\displaystyle \sum _{w=1}^{N_{tw}}P{T}_t\times (\frac{Di{s}_{P-A}}{V_0}}}+{\displaystyle \frac{Di{s}_{P-A}}{V_t^w}})$$

(7)

(3) Assembly stage: The total energy consumption in this stage consists of the energy consumption of the assembly machines, which is composed of two modes of modes: working mode and idle mode. The energy consumption for each mode of the assembly machines is constructed as follows:

Working energy consumption: The energy consumption generated when the assembly machine is in working mode. The working energy consumption in the assembly stage is calculated as follows:

$$E{C}_A^{st3}={\displaystyle \sum _{k=1}^{N_K}{\displaystyle \sum _{a=1}^{N_{Ak}}P{P}_k}}\times P({\pi }_a^k)$$

(8)

Idle energy consumption: This refers to the energy consumption generated between two consecutive operations on the same assembly machine. The idle energy consumption of the machine in the assembly stage is calculated as follows:

$$E{C}_I^{st3}={\displaystyle \sum _{k=1}^{N_K}{\displaystyle \sum _{a=2}^{N_{Ak}}P{I}_k}}\times IP({\pi }_{a-1}^k,{\pi }_a^k)$$

(9)

The total energy consumption in the assembly stage is calculated as follows:

$$E{C}^{st3}=E{C}_A^{st3}+E{C}_I^{st3}$$

(10)

The total energy consumption of the three stages is calculated as follows:

$$TEC(\pi )=E{C}_{st1}+E{C}_{st2}+E{C}_{st3}$$

(11)

Optimization objective: To find an optimal sequence in the set of all three-stage sequences $\prod$($\pi =[{\pi }^{st1},{\pi }^{st2},{\pi }^{st3}]\in \prod$) that simultaneously minimizes the maximum completion time and total energy consumption. The specific details are as follows:

$$C_{\max }(\pi )=\max \{CT({\pi }_a^k)\},\forall a,k$$

(12)

$$TEC(\pi )=E{C}_{st1}+E{C}_{st2}+E{C}_{st3}$$

(13)

$$\{\begin{array}{l}{C}_{\max }({\pi }^{*})=\min \{C_{\max }({\pi }^{st1},{\pi }^{st2},{\pi }^{st3})\},TEC({\pi }^{*})=\min \{TEC({\pi }^{st1},{\pi }^{st2},{\pi }^{st3})\}\\ {\pi }^{*}=\arg \left\{\begin{array}{c}{C}_{\max }({\pi }^{st1},{\pi }^{st2},{\pi }^{st3})\\ TEC({\pi }^{st1},{\pi }^{st2},{\pi }^{st3})\end{array}\right\}\to \min ,\forall ({\pi }^{st1},{\pi }^{st2},{\pi }^{st3})\in \prod \end{array}$$

(14)

3. Encoding and Decoding Design

As described in Section 2.1, the problem characteristics of EEMsMlAJSP result in an exceptionally large solution space. If full encoding is applied to each stage, although the entire solution space can be effectively represented, it also results in the algorithm convergence being slowed and invalid solutions generated during the algorithm iteration process needing to be repaired, thus affecting the search efficiency of the algorithm. Heuristic rules can effectively trim the solution space while maximizing the discovery of near-optimal solutions to the problem, thereby the convergence speed of the algorithm is improved. Therefore, based on the coupling properties between the stages of EEMsMlAJSP [17,18], only the processing operations in the processing stage are encoded, with the assignment of processing machines and the sub-problems in subsequent stages being decoded using heuristic rules.

3.1. Encoding

In this study, the solution to the problem is represented using an operation-based encoding method. The operations of a job are indicated by number of the job, and the number of times that the job number appears represents which operation it is for that job. Due to the design order requirements of the products considered in this paper, and the possibility of identical jobs existing among different products, product information needs to be read sequentially to determine the required jobs for the product, and each job is assigned a unique code for identification. The pseudocode for the encoding is shown in Algorithm 1.

Algorithm 1 Encoding Method

To better illustrate the encoding process, we use an example where the demand for product 1 is 2 and the demand for product 2 is 1. Product 1 consists of components 1 and 2, while product 2 consists of components 3 and 4. Component 1 is composed of job 1, component 2 is composed of jobs 2 and 3, component 3 is composed of jobs 2 and 4, and component 4 is composed of jobs 5, 6 and 7. Each job is processed through corresponding operations. The specific encoding method for this example is shown in Figure 2.

3.2. Decoding Rules

To narrow the solution space while maximizing the discovery of near-optimal solutions to the problem, the heuristic decoding rules for allocation in the processing stage and subsequent transportation and assembly stages are as follows:

(1)

Machine assignment in the processing stage

The earliest completion time (ECT) rule is used to assign operations to the machine that can complete them the earliest. If multiple machines can complete the operation at the same time, the minimum energy consumption first (MECF) rule is applied to select the machine with the lowest energy consumption.

(2)

Jobs transportation sequence and vehicle allocation in the transportation stage

Under the constraint of vehicle load capacity, the first ccompleted first transported (FCFT) rule is used to organize the batch transportation of jobs. This means that jobs are loaded into available vehicles in the order of their completion times, ensuring that no single vehicle is overloaded. Subsequent jobs are transported by following vehicles according to the same rule.

(3)

Product assembly sequence and machine assignment in the assembly stage

When the jobs transported to the assembly shop are sufficient for the complete assembly of the product, the assembly is carried out according to the product assembly process. If multiple products simultaneously meet the assembly conditions, the minimum jobs first (MJF) rule is used to determine the assembly sequence. During product assembly, the machine assignment for assembly operations follows the same rules as in the processing stage to determine the assembly machine for each operation.

3.3. Design of Energy-Efficient Strategy Decoding Trigger Mechanism

To address the issue of substantial idle energy consumption caused by the hierarchical coupling characteristics of EEMsMlAJSP, an energy-efficient strategy based on block-moving can reduce machine idle time by making dispersed operations more compact without changing the maximum completion time, thus reducing the idle energy consumption generated by machines.

However, as the problem scale increases, executing the energy-efficient strategy usually takes a long time and consumes significant computational resources. In fact, not all solutions benefit from this strategy in terms of improving the quality of non-dominated solutions. Therefore, we propose an adaptive trigger mechanism for the energy-efficient strategy decoding, which ensures potential high-quality solutions are not lost while improving the search efficiency. The specific operations are as follows:

Step 1: Apply the decoding rules and the active decoding method [19] to the individual to obtain the scheduling results $C_{\max }(\pi )$ and $TEC(\pi )$, and record the idle time of all machines.

Step 2: According to Equation (15) and Equation (16), obtain the upper bound of energy savings $E{S}_{UB}$ and the lower bound of energy consumption $E_{LB}$ respectively for the scheduling scheme after applying the energy-efficient strategy.

Step 3: Determine the relationship between $(C_{\max }(\pi ),E_{LB})$ and the current non-dominated solutions. If it is not dominated by the current non-dominated solutions or if they are mutually non-dominating, proceed to Step 5, otherwise, proceed to Step 4.

Step 4: Use the objective values $C_{\max }(\pi )$ and $TEC(\pi )$ obtained from the decoding rules as the decoding results of the scheduling scheme.

Step 5: Following the product assembly sequence of the scheduling scheme, check each product in reverse order from the last to the first according to the process sequence and determine whether the block-moving condition is met according to Equation (17). If the condition is met, update the completion time of the operation according to Equation (18). Otherwise, check the preceding operation of this operation.

Step 6: If all operations have been checked, use Equation (19) to determine whether the idle intervals meet the shutdown conditions for profit and loss time, and output the results $C_{\max }(\pi )$ and $TEC(\pi )$ of the scheduling scheme after executing the energy-efficient strategy. Then, update non-dominated solution set. Otherwise, return to Step 5.

$$E{S}_{UB}={\displaystyle \sum _{a=2}^{N_A^k}IP({\pi }_{a-1}^k},{\pi }_a^k)\times P{I}_k$$

(15)

$$E_{LB}=TEC(\pi )-E{S}_{UB}$$

(16)

$$P({\pi }_a^k)\le IP({\pi }_l^k,{\pi }_{l+1}^k)$$

(17)

$$CT({\pi }_a^k)=\{\begin{array}{l}\max (CT{I}_i^k),\forall i,{\pi }_a^k\mathrm{is}\mathrm{the}\mathrm{final}\mathrm{assembly}\mathrm{operation}\mathrm{of}\mathrm{the}\mathrm{product}\\ \min \{C({\pi }_{lasta}^k),\max (CT{I}_i^k)\},\mathrm{otherwise}\end{array}$$

(18)

$$T_{off-on}\le E{C}_{off-on}/\stackrel{—}{EC}$$

(19)

Where $IP({\pi }_l^k,{\pi }_{l+1}^k)$ is the idle time after the operation ${\pi }_a^k$ on machine m, $l=a,a+1,\cdots ,N_A^k-1$. $CT{I}_i^k$ denotes the end time of the insertable idle interval for ${\pi }_a^k$ on machine k, ${\pi }_{lasta}^k$ denotes the immediate subsequent operation of ${\pi }_a^k$. $T_{off-on}$ is the time required for machine turn off and on, $E{C}_{off-on}$ is the energy consumption of for machine turn off and on, and $\overline{EC}$ is the energy consumption per unit time of machine. The flow of energy-efficient strategy decoding trigger mechanism is shown in Figure 3.

4. IICA-QL for EEMsMlAJSP

4.1. Initialization of Empires

In IICA-QL, a country ${\pi }^i$ denotes a solution to the problem, where $i=1,2,\cdots ,Ps$, and Ps denotes the population size. The population $Pop(g)$ in generation g is composed of all countries in generation g, denoted as $Pop(g)=\{{\pi }^{g,1},{\pi }^{g,2},\mathrm{...},{\pi }^{g,Ps}\}$. To obtain a high-quality initial population while maintaining population diversity, two heuristic methods are used to initialize the population. $\lambda \times Ps$ countries, with $\lambda$ set to 0.2, are generated using product aggregation (PA) and job aggregation (JA) rules respectively, while the remaining countries are generated randomly. This approach ensures a high starting point for the search while maintaining search breadth.

Subsequently, according to Equation (20) [20], the cost $c_i^1$ of each country in the initial population are calculated in turn, $i=1,2,\cdots ,Ps$. After sorting $c_i^1$ in non-ascending order, the countries with the highest costs are selected as ${{\rm N}}_{im}$ imperialists for ${{\rm N}}_{im}$ empires, while the other countries in the population are designated as colonies. The number of colonies $N_{col}$ is then calculated according to Equation (21). Next, the normalized power of each imperialist $Im{p}_k$($k=1,2,\cdots ,{{\rm N}}_{im}$) is calculated according to Equation (22), and the number of colonies assigned to each colonial country $N{C}_k$ is determined based on their normalized power. Finally, the corresponding number of colonies is randomly selected and assigned to the corresponding empires, and the selected colonies are removed from the original colony population.

$$c_i^g=\underset{l=1,2,\cdots ,Ps}{\max }\{ran{k}_l\}-ran{k}_i+dis{t}_i/{\displaystyle \sum _{j\in {\theta }_{ran{k}_i}}dis{t}_j}$$

(20)

$$N_{col}=PS-N_{im}$$

(21)

$$Po{w}_k=c_k^g/{\displaystyle \sum _{k=1}^{N_{im}}{c}_k^g}$$

(22)

$$N{C}_k=round(Po{w}_k\times N_{col})$$

(23)

Before the iteration begins, an empty external archive set $\Omega$ is established to store the non-dominated set of the initial population, and it is dynamically updated during the iteration process.

4.2. Assimilation Based on Q-Learning Adaptive Adjustment of Operator Parameter

Assimilation is the process by which colonies within an empire are moved towards the imperialist, with population evolution achieved through the assimilation of colonies by the imperialist. To ensure convergence speed while maintaining population diversity, a dynamic adaptive adjustment strategy for operator parameter based on Q-learning is designed. First, two different operator co-evolution operations are employed. Then, a Q-learning-based adaptive adjustment strategy for the UX operator parameter is implemented to enhance the convergence and diversity of the algorithm. Finally, an elite retention strategy is used to update the individuals of each empire, thereby improving the convergence speed. The pseudocode is shown in Algorithm 2.

Algorithm 2 Assimilation based on Q-learning adaptive adjustment of operator parameter

4.2.1. Adaptive Adjustment of UX Operator Parameter Based on Q-Learning

During the assimilation process, the uniform crossover (UX) operator [21] and the two-point Crossover (TPX) operator [22] are used, respectively applied to the movement of colonies towards the imperialist and the external archive set. Among these, the UX operator is used as the main evolutionary operator to guide the evolutionary process of population and can also be applied for local search [21], playing an important role in the quality of population evolution. The parameter of the UX operator determines the number of high-quality genes passed from parents to offspring. In the early stages of evolution, a larger parameter is helpful for the algorithm converge, while in the later stages of evolution, a larger parameter can result in poor population diversity, causing the algorithm to fall into local optima. Therefore, a Q-learning-based adaptive adjustment strategy for the UX operator parameter is devised to dynamically adjust the evolutionary process of population, enabling the algorithm to converge more efficiently to the Pareto front.

State: The state of the entire population is determined by evaluating the convergence metric $\Delta GD$ of the non-dominated solution set compared to the reference set and the diversity metric $\Delta Ratio$, which represents the coverage ratio of the solution space. Here, $\Delta GD$ represents the generational distance difference of the non-dominated solutions obtained compared to the Pareto front $P{F}^{*}$ between generation t and t+1. $P{F}^{*}$ is a preset reference point, and $\Delta Ratio$ represents the proportion of subregions occupied by the non-dominated solutions in the current generation in the objective space. The objective space is divided into several non-overlapping regions by weight vectors [23]. $\Delta GD$ indicates the degree of convergence of generation t+1 compared to generation t, the smaller the value of $\Delta GD$, the better the convergence. $\Delta Ratio$ represents the change in non-dominated set diversity, with $\Delta Ratio<0$ indicating a decrease in diversity and $\Delta Ratio>0$ indicating an increase in diversity. The calculations for $\Delta GD$ and $\Delta Ratio$ are given by Equation (24) and Equation (25) respectively, where $\left|P{F}_t\right|$ and $\left|{\psi }_t\right|$ represent the number of non-dominated solutions in generation t and the proportion of the objective space they occupy. During the evolutionary, there are four possible states: (1) $\Delta GD>=0,\Delta Ratio>0$;(2)$\Delta GD>=0,\Delta Ratio<=0$;(3)$\Delta GD<0,\Delta Ratio>0$;(4) $\Delta GD<0,\Delta Ratio<=0$.

$$\Delta GD=(\sqrt{{\displaystyle \sum _{x\in P{F}_{t+1}}\underset{y\in P{F}^{*}}{\min }d{(x,y)}^2}}/\left|P{F}_{t+1}\right|)-(\sqrt{{\displaystyle \sum _{x\in P{F}_t}\underset{y\in P{F}^{*}}{\min }d{(x,y)}^2}}/\left|P{F}_t\right|)$$

(24)

$$\Delta Ratio=\left|{\psi }_{t+1}\right|-\left|{\psi }_t\right|$$

(25)

Action: The action set $\mathbb{A}$ is composed of ten actions, representing ten levels of the UX operator parameter. By determining the current state of the population, the action with the highest Q-value in the Q-table is selected based on the ε-greedy strategy.

ε-greedy: In the initial iterations, Q-values usually lack prior knowledge, making it easy for the agent to fall into local optima during learning. Therefore, a ε-greedy strategy based on sigmod is used to guide the learning of agent. The mathematical expression for ε-greedy is as follows:

$$\epsilon =0.1+（\epsilon -0.1）\times {\displaystyle \frac{1}{1+e^{(-0.2\times (g-50))}}}$$

(26)

where g represents the number of iterations.

Reward: By calculating the values of $\Delta GD$ and $\Delta Ratio$ after executing an action, the agent receives a corresponding reward and updates the Q-value of the action in the Q-table for the corresponding state according to Equation (27). If $\Delta GD\ge 0,\Delta Ratio>0$, then reward = 1; if $\Delta GD\ge 0,\Delta Ratio<=0$, then reward = -2; if $\Delta GD<0,\Delta Ratio>0$, then reward = 4; if $\Delta GD<0,\Delta Ratio\le 0$, then reward = 2.

$$Q（s_t,a_t)=Q（s_t,a_t)+\alpha [rewar{d}_{t+1}+\gamma \underset{a\in \mathbb{A}}{\max }Q(s_{t+1},a)-Q（s_t,a_t)]$$

(27)

where $\alpha$ and $\gamma$ are the learning rate and discount factor respectively, and $\mathbb{A}$ is the set of available actions.

Figure 4 provides an example. In generation t, the population is in state S₃ and the action A₂ is chosen. After this action is executed, the population is transitioned to state S₁ in generation t+1. In this state, the action A₁₀ is chosen, and in generation t+2, the population is transitioned to state S₄.

4.2.2. Update Strategy

To achieve rapid convergence of the population, an elite retention strategy is employed to select the colonies that need to be retained after the assimilation process. The specific operations are as follows: individuals generated during the assimilation process of imperialist $Im{p}_k$ are placed into the colony offspring population $Offsprin{g}_k$. After the assimilation operation within the empires is completed, the original colony population is merged with the colony offspring population to form a new population $SubPo{p}_{\mathrm{k}}$. Then, the individuals in $SubPo{p}_{\mathrm{k}}$ are non-dominated sorted and their crowding distances are calculated. A number of $N{C}_{\mathrm{k}}$ colonies are selected as the next generation of colonies. Additionally, for new solutions generated by the assimilation of the imperialist, if they are not dominated by the old solutions, the original imperialist is replaced by the new solutions; otherwise, no replacement is made.

4.3. Revolution Based on Hyper-Heuristic Variable Neighborhood Search

The colonial revolution is a process of further exploring and exploiting the neighborhood of superior colonies after completing the assimilation operation. In this paper, various neighborhood search (VNS) operations are designed, and high-level information constructed from a two-dimensional probability matrix is used to guide the selection of corresponding low-level variable neighborhood search operation sets. The selection of variable neighborhood search operators is dynamically adjusted by the high-level strategy based on the feedback information received from the environment after each application of the low-level variable neighborhood search operation set, thereby improving the search performance of the algorithm.

4.3.1. Neighborhood Search (NS) Operations

To implement an effective perturbation strategy for individuals, two single-operation-based neighborhood operations and five block-based neighborhood operations are designed as neighborhood search operators.

NS1: Single-operation insertion operation for individual encoding. A random operation is selected from the individual encoding and inserted into another position

NS2: Single-operation exchange operation for individual encoding. Two different operations are randomly selected and their positions are swapped.

NS3: Reversal operation for individual encoding. Two positions are randomly selected, and the segment of operations between these positions is reversed.

NS4: Segment exchange operation for individual encoding. Two segments of operations are randomly selected and their positions are exchanged.

NS5: Segment insertion operation for individual encoding. Two segments of operations are randomly selected and one segment is inserted into another position within the encoding.

NS6: Segment shuffling operation for individual encoding. A segment of operations is randomly selected, and operations within segment are randomly shuffled.

NS7: Segment gene dispersal operation for individual encoding. The individual encoding is divided into several segments, and the consecutive repeated genes within each segment are dispersed to different positions.

Figure 5. Schematic Diagram of Neighborhood Structure.

Figure 5. Schematic Diagram of Neighborhood Structure.

4.3.2. Two-Dimensional Probability Model

Define $NS=\{N{S}^{G,k}(1),N{S}^{G,k}(2),\cdots ,N{S}^{G,k}(l)\}$ as the neighborhood operations sequence executed by the k-th individual in generation G, where l represents the length of variable neighborhood search. The two-dimensional probability model $P^G(x,y)$, shown in Equation (28), is used as the high-level strategy to guide the selection of the low-level the variable neighborhood search operations set.

$$P^G(x,y)=\left[\begin{array}{l}{P}_{l\times l}^G(1,1)\cdots P_{l\times l}^G(1,l)\\ ⋮\ddots ⋮\\ P_{l\times l}^G(l,1)\cdots P_{l\times l}^G(l,l)\end{array}\right]$$

(28)

The initialization of the two-dimensional probability model and its update in the G-th generation are shown in Equation (29) and Equation (30) respectively.

$$P^0(x,y)=1/l,x,y=1,2,\cdots ,l$$

(29)

$$P^G(x,y)=（1+reward)\times P^{G-1}(x,y),x=1,2,\cdots ,l,y=NS(1),NS(2),\cdots ,NS(l)$$

(30)

$$reward=\{\begin{array}{l}0.1,\mathrm{new}\mathrm{solution}\mathrm{is}\mathrm{not}\mathrm{dominated}\mathrm{by}\mathrm{old}\mathrm{solution}\\ -0.1,\mathrm{otherwise}\end{array}$$

(31)

After completing the revolution operations, the cost of each country within the empire is calculated, and the country with the highest cost is selected as the new imperialist, while the original imperialist becomes a colony. $P^G(x,y)$ is updated according to Equation (31).

Algorithm 3 Revolution based on hyper-heuristic variable neighborhood search

4.4. Joint Imperialist Invasion

Dividing the population into multiple subpopulations is helpful for individuals to extensively search different areas, and balances the exploration and exploitation capabilities of the algorithm through co-evolution [24]. To address the issues of poor diversity and premature convergence in the basic Imperialist Competitive Algorithm, a joint imperialist competition mechanism is designed to replace the standard imperialist competition operations. Each empire is regarded as a population that collaboratively searches different regions of the solution space. This mechanism ensures that the collaborative evolutionary search is maintained within each empire and that information sharing and interaction are enhanced between empires, thereby improving diversity and the search efficiency of the algorithm.

The specific operations of the joint imperialist invasion are as follows: The total cost of each imperial group is calculated and ranked according to Equation (32), where $Em{p}_k$ represents the set of colonies owned by $Im{p}_k$, and $\xi$ is the proportional coefficient of the total cost of the empire accounted for by the colonies, which is set to 0.1. The imperialists of the top $N_{im}-1$ empires perform POX crossover [4] with random solutions, and the corresponding number of the worst colonies in the worst-ranked empire are replaced by the newly generated solutions. Through joint imperialist competition, the collaborative evolution within each empire can be maintained, while information sharing and interaction between empires can be enhanced, thus improving diversity and increasing the search efficiency of the algorithm.

$$T{C}_k^g=c_i^g+\xi {\displaystyle \sum _{\sigma \in Em{p}_k}{c}_{\sigma }^g}/N{C}_k$$

(32)

Algorithm 4 Joint imperialist invasion

4.5. Algorithm Description

The flowchart of the proposed IICA-QL is shown in Figure 6. The process of IICA-QL is described as follows.

Step 1: Generate the initial population using heuristic rules and random generation methods. Select $N_{im}$ imperialist countries to form the initial empires.

Step 2: Construct the external archive set and update it based on non-dominated solutions.

Step 3: Select UX operator parameter according to the Q-table, then perform assimilation in each empire and update Q-table. Meanwhile, retain colonies based on the elite retention strategy.

Step 4: Execute high-dimensional strategies constructed by the two-dimensional probability model to select low-level neighborhood search operations sequence. Update the probability model based on feedback information.

Step 5: Perform joint imperialist competition to achieve collaborative evolution of empires.

Step 6: If the termination condition is met, stop the search; otherwise, return to Step 3.

5. Experimental Results and Analysis

The proposed IICA-QL is implemented in MATLAB R2021a. The experimental configuration is an Intel Core i7-13700K and CPU @3.40GHz PC with 32GB memory. Since no prior research has been conducted on EEMsMlAJSP, the test data are generated randomly. The problem scale is denoted by O/J/A_M/T/K, where O/J/A represent the number of processing operations, jobs, and assembly operations, respectively. M/T/K represent the total number of processing machines, transport vehicles, and assembly machines, respectively. Detailed information about the product structure, order information, and production parameter data used in the experiments can be found from: https://github.com/captain-dyq/EEMsMlAJSP.git. The experimental time $T_{(i)}$ for each problem scale is determined based on the specific problem scale, where $T_{(i)}=O\times M+J\times T+A\times K$.

5.1. Performance Metrics

To comprehensively evaluate the convergence and diversity of the approximation Pareto fronts set obtained by the algorithm, the following three commonly used metrics are employed to measure the proposed IICA-QL.

(1)

Solution set Coverage (C-metric). This metric measures the dominance relationship between two solution sets, and is defined as follows:

$$\mathrm{C}(A,B)={\displaystyle \frac{|x\in B|\exists y\in A:y\prec x|}{\left|B\right|}}$$

(33)

Where A and B are two non-dominated solution sets obtained by the algorithm. C(A,B) represents the proportion of solutions in set B that are dominated by at least one solution in set A. The larger the value of C(A,B), the better the performance of A.

(2)

Generational Distance (GD).

The proximity of the Pareto front PF obtained by the algorithm to the optimal Pareto front PF* is evaluated, where PF* is the true non-dominated solution set among the non-dominated sets obtained by all algorithms. For each solution x in PF, the nearest solution y in PF* is found, and their Euclidean distance is calculated. GD is the average of these shortest Euclidean distances, which is described as Equation (34). A smaller value of GD indicates better convergence and that the found PF is closer to the optimal Pareto front.

$$GD={\displaystyle \frac{\sqrt{{\displaystyle \sum _{x\in PF}\underset{y\in P{F}^{*}}{\min d{(x,y)}^2}}}}{\left|PF\right|}}$$

(34)

(3)

Hypervolume (HV)

The HV metric is a univariate metric that represents the volume of the region in the objective space enclosed by the non-dominated solution set obtained by the algorithm and the reference point, which is set to (1.01,1.01) [25]. The hypervolume metric is a comprehensive metric for evaluating the convergence and diversity of the non-dominated solution set obtained by the algorithm. The larger the HV value, the better the overall performance of the algorithm.

$$HV(E)=Leb({\cup }_{x\in E}[f_1(x),r_1]\times \mathrm{...}\times [f_m(x),r_m])$$

(35)

Where Leb(·) represents the Lebesgue measure and m represents the number of objectives in the problem.

To eliminate the impact of different dimensions, the solutions in the solution set are normalized according to the following formula [19]:

$$f_i^{\text{'}}={\displaystyle \frac{f_i-f_i^{\min }}{f_i^{\max }-f_i^{\min }}},i=1,2$$

(36)

5.2. Parameter Settings

To determine the optimal values for the parameters in IICA-QL, the Design of Experiment (DOE) method is used for parameter setting. IICA-QL has a total of 7 parameters, each with three levels to choose from. The orthogonal array L27(3⁷) is selected, and the meanings and optional parameter levels for each parameter are shown in Table 2. Experiments are conducted using a medium-sized experimental scale (82/42/71_3/3/3), with each parameter combination run 10 times and the runtime being $T_{(i)}=O\times M+J\times T+A\times K$. The hypervolume (HV) metric is calculated after each run for each parameter combination, and the average of the ten HV values is used as the response value (RV). Additionally, the average RV and influence of each parameter level are shown in Table 3. According to Table 3, the optimal parameter set is $Ps$ = 100， N_im = 3, P_a = 0.8, $P_r$= 0.1, $\alpha$= 0.2, $\gamma$= 0.9, $\epsilon$= 0.6.

5.3. Effectiveness of the Adaptive Trigger Mechanism for the Energy-Efficient Strategy

To verify the effectiveness of the adaptive trigger mechanism for the energy-efficient strategy, IICA-QL (denoted as A) is compared with IICA-QL executing the energy-efficient strategy on Pareto solutions (IICA-QL-PFE, denoted as B) and IICA-QL executing the energy-efficient strategy on all individuals (IICA-QL-AE, denoted as C). Using the same computation time as the termination condition, each algorithm runs independently 10 times on 18 test instances. The average values of C-metric, GD, and HV obtained are shown in Table 4, with the best value for each metric highlighted in bold. The runtime of each algorithm on each scale is recorded in Table 5.

From the C-metric results in Table 4 and Figure 7, it can be seen that in all 18 test instances, C(A,B) is greater than C(B,A), and in 17 instances, C(A,C) is greater than C(C,A). This indicates that the non-dominated solutions obtained by IICA-QL generally dominate those obtained by IICA-QL-PFE and IICA-QL-AE. From the GD and HV metrics, it is evident that in most test instances, the values for IICA-QL are better than those for IICA-QL-PFE and IICA-QL-AE, and the values for IICA-QL-AE are better than those for IICA-QL-PFE. This shows that the adaptive trigger mechanism for the energy-efficient strategy improves the search quality of the algorithm. Additionally, from Table 5 and Figure 8, it can be seen that the CPU time of IICA-QL is comparable to that of IICA-QL-PFE and significantly better than that of IICA-QL-AE, indicating that the adaptive trigger mechanism for the energy-efficient strategy improves the search efficiency of the algorithm. Therefore, the adaptive trigger mechanism for the energy-efficient strategy effectively enhances the performance of the algorithm.

5.4. Effectiveness of Operator Parameter Adaptive Adjustment Based on Q-Learning

To verify the effectiveness of dynamic adaptive adjustment of operator parameter based on Q-learning, IICA-QL (denoted as A) is compared with IICA-QL without this strategy (IICA-nQL, denoted as B), where the UX parameter is set to 0.55. Using the same computation time as the termination condition, each algorithm runs independently 10 times on 18 test instances. The average values of C-metric, GD, and HV obtained are shown in Table 6, with the best value for each metric highlighted in bold.

From the C-metric results in Table 6 and Figure 9, it can be seen that in 17 out of 18 test instances, C(A,B) is greater than C(B,A), and C(A,B) shows an increasing trend. This indicates that the non-dominated solutions obtained by IICA-QL generally dominate those obtained by IICA-nQL, and the advantage becomes more pronounced as the problem scale increases. From the GD and HV metrics, it is evident that in all test instances, the values for IICA-QL are better than those for IICA-nQL, indicating that IICA-QL has better convergence and diversity compared to IICA-nQL. Therefore, the dynamic adaptive adjustment strategy for operator parameter based on Q-learning effectively enhances the performance of the algorithm.

5.5. Comparison with Other Algorithms

Since no existing research has been conducted on solving EEMsMlAJSP, to verify the effectiveness of the algorithm, IICA-QL is compared with algorithms BPBMO [22], PSO-GA[26], and KBOA [17], which solve similar problems to EEMsMlAJSP. The parameters for BPBMO, PSO-GA, and KBOA are chosen as specified in their respective references. Each of the four algorithms is independently run 10 times for the runtime determined by the problem scale, and the average values of the evaluation metrics are calculated. The C-metric, GD, and HV values for each algorithm across 18 problem scales are shown in Table 7 and Table 8, with the best value for each metric highlighted in bold.

As seen in Table 7 and Figure 10, the GD and HV values obtained by IICA-QL are superior to those of BPBMO, PSO-GA, and KBOA in all test instances. Specifically, the performance of IICA-QL in terms of the GD metric indicates that it can find solutions very close to the Pareto optimal front in the solution space, demonstrating a significant advantage in convergence. From the HV metric, it is evident that the values obtained by IICA-QL are much higher than those of the other three algorithms, indicating that IICA-QL not only finds high-quality solutions but also maintains solution diversity, covering a broader Pareto front region. Therefore, the effectiveness of IICA-QL in solving EEMsMlAJSP is validated.

To further illustrate the performance of the compared algorithms, Figure 11 presents the Pareto results obtained by four algorithms on small, medium, and large-scale instances, respectively. As shown in the figure, the Pareto front obtained by IICA-QL outperforms those of the other three algorithms across different instances, indicating that IICA-QL is able to achieve non-dominated solutions with lower $C_{\max }$ and TEC compared to BPBMO, PSO-GA, and KBOA.

6. Conclusions

This study has addressed the multi-stage, multi-level assembly job shop scheduling problem, aiming to minimize the maximum completion time and total energy consumption. An improved imperialist competitive algorithm based on Q-learning (IICA-QL) has been proposed. First, the EEMsMlAJSP problem model has been established. Second, encoding and decoding strategies based on the problem's characteristics have been designed, and an energy-saving decoding strategy trigger mechanism has been proposed to improve solution quality and search efficiency. Then, a dynamic adaptive adjustment strategy for the assimilation operator parameter has been achieved through Q-learning, enhancing convergence speed while ensuring population diversity. Finally, Experimental results have demonstrated the effectiveness of the proposed algorithm.

Future research will further consider the following aspects: (1) incorporating more complex constraints and optimization objectives widely present in actual production to make the problem model more aligned with real production needs; (2) applying other algorithms to solve EEMsMlAJSP, combining optimization algorithms with machine learning techniques to enhance algorithm performance; (3) improving the diversity of the algorithm's search to provide a better-distributed non-dominated solution set.