1 Introduction

Aggregate production planning (APP) is an operational activity that provides an aggregate plan for production processes. The aim of APP is to set overall output levels to face with fluctuating demands. APP is a decision-making process to determine the best way to utilize resources to meet forecasted demand. Since the introduction of APP problem in 1950s, it has been studied vastly by many researchers. The interest in APP has a root in the ability that it provides to companies for effective control of production and inventory costs as the two substantial portions of the overall cost of manufacturers [1].

Aggressive industrial market competition makes enterprises must face the changing market environment. The diversity of products and fierce competition make the stability of manufacturing industry and supply chain more important than ever before. So, enterprises should be to arrange the production plan reasonably to improve the enterprise's adaptability and reduce the production cost. Many researchers have studied APP problem vastly, because it can help to determine the problems of raw material procurement plan, the number of workers, overtime and inventory levels, etc. [1,2,3].

According to the number of objective functions in the models proposed in the literature, APP models can be classified into two categories: single-objective APP problem and multi-objective APP problem.

1.1 Single-Objective Function APP

Over the last decades, numerous single-objective APP models have been studied. Techawiboonwong and Yenradee [4] presented a mathematical model for a multiple product-types APP to minimize the comprehensive cost, and the spreadsheet-solver technique was used as a tool to solve the model. Tang et al. [5] studied a multi-product APP problem mainly considering the production capacity, capital level and inventory capacity of the enterprise. Paiva and Morabito [6] presented a mixed-integer programming model to maximize the total variable revenue in the APP of sugar mills considering the inventory balance, the capacity constraint, the small bucket constraint, the compatibility constraint, etc. Zhu et al. [7] established an optimization model of a multi-period and multi-produce APP with the production capacity and demands constrains to minimize the total cost. Zhai et al. [8] established an optimization model of mass customization-oriented remanufacturing APP to maximize the profit considering the supply and demand, the inventory capacity, the production capacity and transport capacity. Ramezanian et al. [9] developed a mixed-integer linear programming (MILP) model for general two-phase APP systems, and implemented a genetic algorithm and tabu search for solving this problem. Erfanian and Pirayesh [10] introduced a MILP model for the integration of APP to minimize the total production and maintenance planning costs. Hossain et al. [11] developed an APP model to minimize the total costs of inventory levels, overtime, labor levels, and capacity, etc., and genetic algorithm optimization (GAO) approach and Big M method were used for solving a real time multi-product, multi-period APP decision problem. Jang and Chung [12] proposed a robust optimization approach to the APP problem under implementation errors related to hiring and layoff.

The single objective of most studies used for APP model is to minimize the total cost or maximize the total revenue of the system. In the actual production process, the production planning management of enterprises not only needs to consider from the cost, revenue, quality and other single aspects, but also needs to take into account customer satisfaction, the overall cost of supply chain, staff stability and other aspects.

1.2 Multi-objective Function APP

In addition to the total cost or the revenue, there are other objectives always considered, such as total revenue, service level, the changing of the workforce level, and customer satisfaction subject to bounds on inventory, backorder, subcontracting.

Wang and Liang [13] proposed a multi-objective APP model including three objective functions of the total costs, the carrying and back ordering costs, and the changing workforce level. Leung and Chan [14] developed a multi-objective APP model to maximize the profit, minimize the repairing costs, and maximize machine utilization. Alehashem et al. [15] addressed a multi-objective mixed-integer nonlinear programming model considering two conflicting objectives. Sadeghi et al. [16] implemented a goal programming approach to minimize the total costs, carrying and back ordering costs, and the rate of changes in workforce level. Madadi and Wong [3] developed a multi-objective integer linear programming APP model to minimize the cost and maximize the customer service level. Modarres and Izadpanahi [17] proposed a linear multi-objective optimization APP model to minimize operational costs, energy costs and carbon emission, and a robust optimization approach is applied. Jamalnia et al. [18] proposed a novel decision model to APP decision-making problem considering total revenue, total production costs, total labor productivity costs, optimum utilization of production resources and capacity and customer satisfaction. Mehdizadeh et al. [19] developed a bi-objective optimization model for an APP problem to maximize the profit by improving learning and reducing the failure cost of the system and minimize the costs associated with repairs and deterioration. Djordjevic et al. [20] proposed a fuzzy model for optimal APP which minimizes the total time required to manufacture products, store them in a warehouse and prepare for delivery to customers. Rasmi et al. [21] presented a multi-objective APP model to analyze economic, social, environmental, and cultural pillars inclusively.

Since APP problem always involves several objectives, multi-objective programming has been widely used in this area. A review of the literature on APP reveals that most of the existing research on APP mainly considers the constraints on the balance equation for production, inventory capacity, inventory and demand, production capacity and labor capacity. Almost all researchers have considered the total production cost or sales profit as the first objective to be minimized, and minimization of inventory and backorder level are the other main objectives that have been taken into consideration. The main methods applied to deterministic management science techniques for APP problem are: linear programming, piecewise linear programming, nonlinear/quadratic programming, etc.

Because of the conditions of labor market, the proportion of labor cost is increasing. Decision makers pay great attention to the labor cost. Furthermore, the fluctuation of demand will lead to the change of enterprise employment, and the change of workers will cause the change of labor intensive and labor cost. So, the stability of employment is very important to APP problem. With the development of manufacturing industry from large-scale production of a single product to intensive production mode of multi-product and multi-stage, more dynamic characteristics of production need to be considered in practical application. Therefore, the APP model should include multiple factors related to production and inventory, such as production costs, labor productivity costs, overtime production costs, inventory costs and the changing of the workforce level. In order to improve the adaptability of enterprises, especially the stability of enterprise production and employment, and reduce the production cost, it is necessary for the enterprise to arrange the production plan comprehensively considering various factors.

Due to the large scale and many factors of APP problems, it is very important to an effective algorithm. Meta-heuristics algorithm has been proved to be an effective solution algorithm. Among the meta-heuristics, genetic algorithms (GAs) [9, 11, 19, 22,23,24], particle swarm optimization (PSO) [12, 25, 26], tabu search (TS) [9, 27, 28] and harmony search algorithm (HSA) [19, 29, 30] have been used to deal with APP models [31]. One of the widely used algorithms among meta-heuristics proposed to solve multi-objective APP problems is the Pareto-based non-dominated sorting genetic algorithm called NSGA-II (non-dominated sorting genetic algorithm-II).

This paper develops a bi-objective APP model of multi-product, multi-stage to minimize the total production costs and instability in the work force, considering various factors, the relationship of raw materials, inventory cost and product demand. As the complexity involved in this APP model, a local search-based GA (LS-GA) is designed to improve the efficiency of the algorithm. GA has the characteristics of random multi-point search and implicit parallelism, and it is not easy to fall into the local optimal solution. It has become a common method to solve complex combinatorial optimization problems, while a simple genetic algorithm is prone to premature and has a slow convergence. Local search (LS) algorithm can improve the search ability in the solution space. Due to the combination and large scale of APP problem, a GA is designed to solve the model combined with the local search algorithm based on the minimum-cost flow.

The remainder of this paper is organized as follows: Some mathematical notations are defined in Sect. 2.1. The constraints and objective function are discussed in Sects. 2.2 and 2.3, respectively, and then the bi-objective APP problem model is established in Sect. 2.4. A GA based on the local search algorithm is proposed in Sect. 3. The proposed algorithm is tested via some numerical examples in Sect. 4. Finally, Sect. 5 provides concluding remarks and suggests some future research directions.

2 Bi-objective Model for APP

To adapt to the rapid changes in market demand, the mode of multi-product small-batch and short-period production is becoming more widely applied in the manufacturing industry, and in the meanwhile, the satisfactory of customer and production economy are attached more attention. For most manufacturing enterprises, production planning can be divided into several production planning stages, including raw material procurement, product production, inventory control and labor arrangement. Each stage will change in different periods, such as seasonal demand, price fluctuation of raw materials, etc. Furthermore, the fluctuation of demand will lead to the change of enterprise employment, and the change of workers will cause the change of labor intensity and labor cost. The main objective of the paper is to find out the optimization scheme of labor arrangement, production plan, inventory and raw material procurement in the planning period by considering the relationship between raw materials, stability of worker, inventory cost and product demand under the condition of deterministic demand.

So, a bi-objective optimization model is developed in this paper for APP problem considering the relationship of raw materials, inventory cost and product demand. The first objective function minimizes total production costs. The second objective function is designed to increase the stability in the work force. These two objective functions are formulated in Sect. 2.3 and the constraints of the supply and demand relationship, inventory capacity, production capacity, labor capacity, etc. are analyzed in Sect. 2.2 based on the mathematical notations defined in Sect. 2.1.

The assumptions about the APP problem are as follows:

  1. 1.

    The raw materials purchased in each period are sufficient and no surplus.

  2. 2.

    The demand for the products in each period is known constant and the delivery of products cannot be delayed.

  3. 3.

    All workers are competent for the production of various products, and the skills of workers have no difference.

  4. 4.

    The inventory cost of one product is same in each period.

2.1 Notations

Let \(t = 1,2, \ldots ,T\) be the production planning period set, where T is the number of periods in planning horizon. Let \(i = 1,2, \ldots ,I\) be the product category set, where \(I\) is the number of product category, and \(j = 1,2, \ldots ,J\) be the raw material category set, where \(J\) is the number of raw material category.

2.1.1 Notations of product

\({\text{PD}}_{it}\): Demand for product i in period t (units);

\({\text{PP}}_{it}\): Production quantity of product i in period t (units);

\({\text{PC}}_{i}\): Unit production cost of product i;

\(C_{1t}\): Total production cost in period t;

\({\text{PR}}_{it}\): The unit cost of raw material for product i in period t;

\(C_{2t}\): Total raw material cost in period t;

\({\text{PW}}_{i}\): Labor hours required for producing a unit product i;

\({\text{PN}}_{it}\): The production capacity of equipment for product i in period t.

2.1.2 Notations of Inventory

\({\text{CI}}_{it}\): The inventory of product i in period t (units);

\({\text{CK}}_{i}\): Inventory cost of product i;

\({\text{CN}}_{i}\): The inventory capacity of product i;

\(C_{3t}\): Total inventory cost in period t.

2.1.3 Notations of Raw Material

\(R_{ij}\): demand for raw material j for producing unit product i;

\({\text{RM}}_{jt}\): Total demand for raw material j in period t;

\({\text{RC}}_{jt}\): The price of raw material j in period t.

2.1.4 Notations of Worker

\(W_{t}\): The number of workers in period t;

\({\text{WH}}_{t}\): The number of workers hired in period t;

\({\text{WL}}_{t}\): The number of workers laid off in period t;

\({\text{WHC}}\): Training cost to hire one worker;

\({\text{WR}}\): Maximum labor hours for regular time in each period;

\({\text{WO}}\): Maximum labor hours for overtime in each period;

\({\text{WRT}}_{{\text{t}}}\): Total labor hours for regular time in period t;

\({\text{WOT}}_{t}\): Total labor hours for overtime in period t;

\({\text{WRC}}\): Regular time labor cost per hour;

\({\text{WOC}}\): Overtime labor cost per hour;

\(C_{4t}\) Total labor cost in period t.

2.2 Analysis of Constraints

Production and inventory balance:

$$ {\text{CI}}_{it} = {\text{CI}}_{it - 1} + {\text{PP}}_{it - 1} - {\text{PD}}_{it - 1} . $$
(1)

Workforce balance:

$$ W_{t} = W_{t - 1} + {\text{WH}}_{t} - {\text{WL}}_{t} . $$
(2)

Constraints (3) and (4) represent the inventory capacity and production capacity in each time period.

$$ 0 \le {\text{CI}}_{it} \le {\text{CN}}_{i} , $$
(3)
$$ 0 \le {\text{PP}}_{it} \le {\text{PN}}_{it} . $$
(4)

Constraint (5) ensures that the initial inventory level and volume products in each period should be equal or greater than the market demand.

$$ {\text{PP}}_{it} + {\text{CI}}_{it} \ge {\text{PD}}_{it} . $$
(5)

Constraint (6) ensures that the total production time of products cannot exceed the total time of workers.

$$ \sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } \le W_{t} ({\text{WR}} + {\text{WO}}). $$
(6)

2.3 The Objective Function

The total production cost and stability in the workforce are mainly considered in this paper, where the total production cost mainly includes the product production cost, raw material cost, product inventory cost and labor cost.

The total production cost in period t can be calculated by

$$ C_{1t} { = }\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} } \cdot {\text{PC}}_{i} . $$
(7)

Total demand for raw material j in period t can be calculated by

$$ {\text{RM}}_{jt} = R_{ji} \cdot {\text{PP}}_{it} . $$
(8)

The unit cost of raw material for product i in period t can be expressed as follows:

$$ {\text{PR}}_{it} { = }\sum\limits_{j = 1}^{J} {R_{ij} \cdot {\text{RC}}_{jt} } . $$
(9)

Therefore, the total raw material cost in period t can be written as:

$$ C_{2t} { = }\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PR}}_{it} } { = }\sum\limits_{i = 1}^{I} {\sum\limits_{j = 1}^{J} {{\text{PP}}_{it} \cdot R_{ij} \cdot {\text{RC}}_{jt} } } . $$
(10)

The total inventory cost in period t is

$$ C_{3t} { = }\sum\limits_{i = 1}^{I} {{\text{CI}}_{it} \cdot {\text{CK}}_{i} } . $$
(11)

If the working hours required for production are less than the maximum labor hours for regular time of workers, \(\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } \le W_{t} \cdot {\text{WR}}\), the total labor hours for regular time in period t would be \(\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} }\) and the total labor hours of overtime would be 0. Otherwise, the total labor hours for regular time in period t would be \(W_{t} \cdot {\text{WR}}\) and the total labor hours of overtime would be \(\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } - W_{t} \cdot {\text{WR}}\). Thus, the total labor hours of regular time can be written as:

$$ {\text{WRT}}_{t} = \min \left\{ {\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } ,\;W_{t} \cdot {\text{WR}}} \right\}. $$
(12)

and the total labor hours for overtime can be expressed as:

$$ {\text{WOT}}_{t} = \max \left\{ {\;\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } - W_{t} \cdot {\text{WR}},0} \right\}. $$
(13)

The total labor cost in period t is, thus, calculated by

$$ C_{4t} {\text{ = WH}}_{t} \cdot {\text{WHC}} + W_{t} \cdot {\text{WI}} + {\text{WRT}}_{t} \cdot {\text{WRC}} + {\text{WOT}}_{t} \cdot {\text{WOC}}{.} $$
(14)

So, the first objective function to minimize the total production cost is as follows:

$$ Z_{1} = \sum\limits_{t = 1}^{T} {\left( {C1_{t} + C2_{t} + C3_{t} + C4_{t} } \right)} . $$
(15)

As the main factor of enterprise production capacity, the stability of the enterprise's production capacity will be affected by the hiring and layoff of workers. Therefore, it is necessary to minimize the change of workers when making production plans. Hence, in the second objective function, the sum of changes in the number of workers is used to measure its stability.

$$ Z_{2} = \sum\limits_{t = 1}^{T} {\left( {{\text{WH}}_{t} + {\text{WL}}_{t} } \right)} . $$
(16)

2.4 The Bi-objective Optimization Model

According to the above analysis, considering the production capacity, inventory capacity, supply and demand relationship and labor level, the bi-objective APP problem model can be formulated as follows:

$$ \begin{gathered} \min \;Z_{1} = \sum\limits_{t = 1}^{T} {\left( {C1_{t} + C2_{t} + C3_{t} + C4_{t} } \right)} \hfill \\ \min \;Z_{2} = \sum\limits_{t = 1}^{T} {\left( {{\text{WH}}_{t} + {\text{WL}}_{t} } \right)} \hfill \\ {\text{s.t.}} \hfill \\ \begin{array}{*{20}l} {{\text{CI}}_{it} = {\text{CI}}_{it - 1} + {\text{PP}}_{it - 1} - {\text{PD}}_{it - 1} ,\;\forall i,t} \hfill \\ {W_{t} = W_{t - 1} + {\text{WH}}_{t} - {\text{WL}}_{t} ,\;\forall t} \hfill \\ {0 \le {\text{CI}}_{it} \le {\text{CN}}_{i} ,\;\forall i,t} \hfill \\ {0 \le {\text{PP}}_{it} \le {\text{PN}}_{it} ,\;\forall i,t} \hfill \\ {{\text{PP}}_{it} + {\text{CI}}_{it} \ge {\text{PD}}_{it} ,\;\forall i,t} \hfill \\ {\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } \le W_{t} ({\text{WR}} + {\text{WO}}),\;\forall t} \hfill \\ {W_{t} ,{\text{PP}}_{it} ,{\text{CI}}_{it} \in N,\;\forall i,t.} \hfill \\ \end{array} \hfill \\ \end{gathered} $$

3 Genetic Algorithm

The genetic algorithms (GA), as powerful and broadly applicable stochastic search and optimization techniques, are widely known types of evolutionary computation methods today. In general, a GA has five basic components: (1) a genetic representation of potential solutions to the problem; (2) a way to create a population (an initial set of potential solutions); (3) an evaluation function rating solutions in terms of their fitness; (4) Genetic operators that alter the genetic composition of offspring (selection, crossover, mutation, etc.); and (5) parameter values that genetic algorithm uses (population size, probabilities of applying genetic operators, etc.) [33]. A simple template for the operation of a GA is summarized as Procedure 1.

figure a

In this section, a GA based on local search is designed to solute the bi-objective APP problem model. The framework of GA mainly includes the following four parts: The first is the representation and the decoding. The second part is the generating of an initial solution. The third part serves as the role of the neighboring by crossover, mutation and local search. The last part is an iterate process by generating new solutions and selecting mechanism. In this section, the representation, decoding and initialization are introduced in Sects. 3.1–3.2, respectively. In Sect. 3.3, the crossover operator, mutation operator and local search based on augmenting cycle algorithm are proposed and, finally, three evaluation methods are analyzed in Sect. 3.4. The flow chart of the proposed LS-GA is shown in Fig. 1.

Fig. 1
figure 1

The flow chart of the proposed LS-GA

3.1 Chromosome Structure

In this study, the production quantity of product for each period \(PP_{it}\) and the number of workers \(W_{t}\), the decision variables, both are positive integers, which are conducive to the implementation of chromosome encoding. So, \({\text{CH}} = (P,W)\) is designed as a chromosome. For sub-chromosome P, an integer I × T matrix is adopted to represent the chromosome of production quantity of product, where I represents the number of product category and T is the number of periods. A gene in this chromosome represents the production quantity (\({\text{PP}}_{it}\)) of product i in period t. For sub-chromosome W, a vector with t elements is adopted to represent the sub-chromosome of \(W_{t}\) for each period.

Figure 2 shows a chromosome with T periods and I product categories, in which gene \(PP_{it}\) represents the production quantity of product and gene \(W_{t}\) represents the number of workers in each period.

Fig. 2
figure 2

A simple example of a chromosome

3.2 Initialization

The initialization is an important process in GA, serving the role of randomly initializing the solution in the feasible region [32]. In what follows, the feasible regions of two decision variables in each period are analyzed.

3.2.1 Analysis of the Production Quantity of Product

The minimum production quantity of each product should meet the demand for the product in the planning period according to the Constraint (5). Since the inventory at the beginning of each period is known, the minimum production quantity can be obtained using the following equation:

$$ {\text{Min}}\;{\text{PP}}_{it} = \max \{ 0,{\text{PD}}_{it} - {\text{CI}}_{it} \} . $$
(17)

Moreover, the maximum production quantity of each product in the planning period should not exceed the production capacity and inventory capacity according to Constraints (3) and (4).

$$ {\text{Max}}\;{\text{PP}}_{it} = \min \left\{ {{\text{CN}}_{i} + {\text{PD}}_{it} - {\text{CI}}_{it} ,{\text{PN}}_{it} } \right\}. $$
(18)

3.2.2 Analysis of the Number of Workers

The total labor hours in each period is \(\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} }\), so the minimum number of workers can be determined as

$$ {\text{Min}}W_{t} = \left\lceil {\left( {\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } } \right)/({\text{WR}} + {\text{WO}})} \right\rceil . $$
(19)

The worker’s number of the current period can be calculated according to the production quantity, since the number of workers in the previous period is known. Thus, the worker’s maximum number can be obtained according to the labor cost and stability in the workforce.

$$ {\text{Max}}\;W_{t} = \max \left\{ {W_{t - 1} ,\left\lceil {\left( {\sum\limits_{i = 1}^{I} {{\text{PP}}_{it} \cdot {\text{PW}}_{i} } } \right)/{\text{WR}}} \right\rceil } \right\}. $$
(20)

Obviously, the obtained chromosomes satisfy the Constraints (3)–(6) in the process of initialization, when the sub-chromosome P and W are produced in the feasible regions of the production quantity and number of workers.

The following algorithm shows the process of initializing a chromosome.

  • Step1. Parameter setting: set T, I, J, \({\text{PD}}_{it}\), \({\text{PW}}_{i}\), \({\text{PN}}_{it}\), \({\text{CK}}_{i}\), \({\text{CN}}_{i}\), \({\text{RC}}_{jt}\), \({\text{WHC}}\), \({\text{WR}}\), \({\text{WO}}\), \({\text{WI}}\), \({\text{WRC}}\), \({\text{WOC}}\).

  • Step2. For t = 1 to T, repeat Step3 to Step4.

  • Step3. For i = 1 to I, calculate the minimum production quantity \({\text{MinPP}}_{it}\) and maximum production quantity \({\text{MaxPP}}_{it}\), and generate a random integer number \(PP_{it} \in [{\text{MinPP}}_{it} {\text{,MaxPP}}_{it} ]\) as a gene of sub-chromosome P. Calculate the inventory of products for next planning period \(CI_{it + 1} { = }CI_{it} + PP_{it} - PD_{it}\).

  • Step4. Calculate the minimum and maximum number of workers \({\text{MinW}}_{t}\),\({\text{MaxW}}_{t}\), and generate a random integer number \(W_{t} \in [{\text{MinW}}_{t} {\text{,MaxW}}_{t} ]\) as a gene of sub-chromosome W.

3.3 Genetic Operators

3.3.1 Crossover Operator

According to the crossing probability \(P_{c} \in (0,1)\), some chromosomes are selected as parents and operated according to the following crossover operators.

3.3.1.1 Partheno Crossover Operator

This crossover operator is only utilizable for sub-chromosome \(P\). Randomly generate two unequal periods \(t_{1} ,t_{2} \in [1,T]\) of the same product for a selected parent chromosome, and then exchange the production quantity. The partheno crossover process is illustrated in Fig. 3.

Fig. 3
figure 3

Example for partheno crossover operator

Not always is a partheno crossover operator able to generate a feasible solution. In this case, a repair procedure is necessary to obtain feasible solutions. These repair methods can be divided into two subcategories: repairing the unfeasible gene of sub-chromosome P and repairing the unfeasible gene of sub-chromosome W. A method is presented below for each subcategory.

  • Repairing unfeasible gene of P: Check the exchanged gene whether in the feasible regions of the production quantity. As shown in Fig. 3, we need to check \({\text{MinPP}}_{21} \le 40 \le {\text{MaxPP}}_{21}\) and \({\text{MinPP}}_{2T - 1} \le 45 \le {\text{MaxPP}}_{2T - 1}\). For unsatisfied genes, regenerate the genes in the feasible regions. That is, generate a random integer number \({\text{PP}}_{it} \in [{\text{MinPP}}_{it} ,{\text{MaxPP}}_{it} ]\) as a gene of P’.

  • Repairing unfeasible gene of W: Check the gene of worker number \(W_{{t_{1} }} ,W_{{t_{2} }}\) whether in the feasible regions of the number of workers. The unsatisfied genes will be regenerated in the feasible regions. That is, generate a random integer number \(W_{t} \in [{\text{Min}}W_{t} ,{\text{Max}}W_{t} ]\).

3.3.1.2 Arithmetic Crossover

Arithmetic crossover is shown in Fig. 4. The method of this operator is as following:

$$ \begin{array}{*{20}c} \begin{gathered} p^{\prime}_{it} = \left\lceil {\lambda p1_{it} + (1 - \lambda )p2_{it} } \right\rceil , \hfill \\ w^{\prime}_{it} = \left\lceil {\lambda w1_{t} + (1 - \lambda )w2_{t} } \right\rceil , \hfill \\ \end{gathered} & {\lambda \in (0,1)} \\ \end{array} , $$
(21)
Fig. 4
figure 4

Illustration of the arithmetic crossover

where \(p1_{it}\),\(p2_{it}\),\(w1_{t}\),\(w2_{t}\) are the genes of the production quantity and worker numbers for selected parents P1 and P2, and \(p^{\prime}_{it}\),\(w^{\prime}_{it}\) are the genes of offspring.

3.3.2 Mutation Operator

A parameter \(P_{m} \in (0,1)\) is defined as the probability of mutation. Some chromosomes are selected to operate the following mutation operators. Two mutation operators are applied to sub-chromosome P and W, respectively.

3.3.2.1 The Production Mutation

This mutation operator is utilizable for sub-chromosome P. Randomly select a gene from P and regenerate this gene in the feasible region. The process of production mutation is shown in Fig. 5.

Fig. 5
figure 5

Illustration of the production mutation

Similarly, the feasibility of genes that selected and subsequent genes need to be examined. Repair methods for the unfeasible gene are the same as the partheno crossover operator.

3.3.2.2 The Mutation for Number of Workers

For sub-chromosome W, randomly select a gene from W and regenerate this gene in the feasible region as shown in Fig. 6.

Fig. 6
figure 6

Illustration of the mutation for number of workers

3.3.3 Local Search Algorithm

Local search algorithm can effectively improve the quality of the current solution [33]. A method is presented below for each sub-chromosome.

3.3.3.1 Local Search Based on Augmenting Cycle Algorithm

Under the condition that the sub-chromosome W remains unchanged, the local search algorithm searches the neighborhood from an initial solution of production quantity to only improve the first objective function \(Z_{1}\).

When the sub-chromosome W is determined, the second objective and the Constraint (2) do not need to be considered. In this case, the problem of optimization of production quantity is a special case of the minimum-cost flow (MCF) problem. Since the total labor cost \(C_{4t}\) is a piecewise function, this problem cannot be solved by network simple method directly. The core idea of solving the MCF problem is augmenting flow in the negative cycle. Negative cycle in cost network is the cycle that the sum of costs of all the edges in the cycle is negative. Augmenting flow in the negative cycle can decrease the cost and improve the quality of solution effectively.

First, we can deal with the production planning as a network model, as shown in Fig. 7. Based on the above idea, if we can find a negative-cost cycle on the network, the cost can be reduced. Let us augment \(\theta_{I}\) units of flow along the cycle in the direction of its orientation. Increases flow on forward arcs by \(\theta_{I}\) units and decreases the flow on backward arcs by \(\theta_{I}\) units, as shown in Fig. 8.

Fig. 7
figure 7

The network model of production planning

Fig. 8
figure 8

Illustration of augmenting flow in the negative cycle

For sub-chromosome P, randomly select two unequal periods \(t_{1} ,t_{2} \in [1,T]\), and a cycle is formed. If the cycle cost is negative, we send \(\theta_{I}\) flow around the cycle in the direction of traversal. If it is positive, we send \(\theta_{I}\) flow around the cycle in the opposite direction. Since the augmenting flow is a multi-product flow, it is necessary to determine the adjustment amount of each product. If the formed cycle is anticlockwise, the maximum adjustment amount of product i can be obtained using the following equation:

$$ {\text{FP}}_{{\text{i}}} = \min \left\{ {{\text{PN}}_{{it_{1} }} - {\text{PP}}_{{it_{1} }} ,{\text{PP}}_{{it_{2} }} ,{\text{CN}}_{i} - {\text{CI}}_{it} |t_{1} < t \le t_{2} } \right\}. $$
(22)

If the formed cycle is clockwise, the maximum adjustment amount can be obtained by using the following equation:

$$ {\text{FP}}_{{\text{i}}} = \min \left\{ {{\text{PN}}_{{it_{2} }} - {\text{PP}}_{{it_{2} }} ,{\text{PP}}_{{it_{1} }} ,{\text{CI}}_{it} |t_{1} < t \le t_{2} } \right\}. $$
(23)

The process of local search based on augmenting cycle algorithm is stated as follows:

  • Step1. Randomly generate two unequal periods \(t_{1} ,t_{2} \in [1,T]\).

  • Step2. Determine the negative-cost cycle formed from two periods \(t_{1} ,t_{2}\).

  • Step3. Calculate the maximum adjustment amount of product i by using the Eq. (22) or (23).

  • Step4. Perform the following operations m times.

    • Step4.1. For each \(i \in I\), generate a random adjustment amount \(\theta_{i} \in [0,FP_{i} ]\).

    • Step4.2. Compute the improving of first objective function \(\Delta Z_{1}\).

    • Step4.3. If \(\Delta Z_{1} < 0\), then adjust the values of corresponding genes.

3.3.3.2 Neighborhood Search for Number of Workers

For sub-chromosome W, randomly select a gene w from W and a neighborhood parameter \(\delta\) to be chosen appropriately, where\(w - \delta ,w + \delta \in {\text{[Min}}W,{\text{Max}}W]\). The iterative improvement algorithm is carried out in the range of \([w - \delta ,w + \delta ]\).

3.4 Fitness Assignment and Selection

The evolutionary algorithm for multi-objective iterates by generating new solutions and selecting mechanism. Therefore, how to evaluate the fitness of individuals and make a select to guide the search to the Pareto optimal set is a very important process [34, 35]. The general evaluation methods of multi-objective individuals include: Non-Pareto methods (such as, \(\varepsilon\)-Constraint Method, Synthesized Objective Value Method, Vector Evaluated Genetic Algorithm), Pareto-based methods (such as NSGA-II, Strength Pareto Evolutionary Algorithm, Pareto Archived Evolution Strategy) and multi-objective Evolutionary Algorithm based on Decomposition (MOEA/D) (such as Multi-Population Genetic Algorithm, Sub-population Evolutionary Algorithm).

To analyze the efficiency of the algorithms, three evaluation methods are employed.

3.4.1 An Entropy-Based Evaluation Method (EBEGA)

Information entropy (IE) is one of the commonly used objective weighting methods [36]. The weights of two objectives can be determined based on the use of the IE method and the synthesized objective value of each individual can be obtained. Furthermore, the selection is made by a roulette-wheel method.

3.4.2 Multi-population Genetic Algorithm (MPGA)

The population of each generation is divided into two subgroups according to the objective functions. Moreover, genetic operators are operated independently on two subgroups. For each subgroup, the selection is done by a roulette-wheel method according to the single objective [37, 38].

3.4.3 NSGA-II

NSGA-II is one of the widely used algorithms among the meta-heuristics proposed to solve multi-objective problems. In this paper, non-dominated sorting and crowding distance-based fitness assignment strategies of NSGA-II algorithm are used to implement the competitive selection [39].

4 Numerical Experiments

Some experimental studies were conducted to compare the performance of the proposed evaluation of NSGA-II with the ones of EBEGA and MPGA.

4.1 Performance Measures

It is well known that multi-objective optimization problems need multiple, uniformly distributed solutions to form a Pareto frontier. We adopt some multi-objective performance measures to evaluate the performances of NSGA-II, EBEGA and MPGA.

  1. 1.

    The average objectives value, \({\text{avg}}( \cdot )\). The smaller the average value, the better the solution quality of the algorithm [40].

  2. 2.

    The number of non-dominated set, \(M_{1}\). The multiple non-dominated solutions are useful for decision makers.

  3. 3.

    Runtime. The runtime of the algorithms is another criterion. Here, average runtime of each algorithm is used.

  4. 4.

    Set coverage measure, \(M_{2}\). \(M_{2}\) was proposed by Miandoabchi, Farahani and Szeto [41] to compare the performance of two algorithms based on the concept of set coverage measure given in Zitzler, Deb and Thiele [42].

    Let \(X^{\prime } ,X^{\prime \prime } \in X\) be two non-dominated sets of different algorithms. The function C mapping the ordered pair \((X^{\prime } ,X^{\prime \prime } )\) into the interval [0, 1] is defined as [42]

    $$ C\left( {X^{\prime } ,X^{\prime \prime } } \right) = \frac{{|\left\{ {x^{\prime \prime } \in X^{\prime \prime } |\exists x^{\prime } \in X^{\prime } ,x^{\prime } \prec x^{\prime \prime } } \right\}|}}{{|X^{\prime \prime } |}}. $$
    (24)

    Note that the dominance rate, \(C\left( {X^{\prime } ,X^{\prime \prime } } \right)\), is not necessarily equal to \(1 - C\left( {X^{\prime \prime } ,X^{\prime } } \right)\). So, based on this set coverage concept, the following equation can be used to express the effect of one algorithm with relative to another algorithm by using the following equation:

    $$ M_{2} = C\left( {X^{\prime } ,X^{\prime \prime } } \right) - C\left( {X^{\prime \prime } ,X^{\prime } } \right). $$
    (25)
  5. 5.

    Mean ideal distance (MID). MID was defined by Zitzler and Thiele [43].

    $$ {\text{MID}} = \sum\limits_{i = 1}^{n} {c_{i} /n} , $$
    (26)

    where n is the number of non-dominated sets and \(c_{i} { = }\sqrt {(Z_{1} (i)/5000 - o_{1} )^{2} + (Z_{2} (i) - o_{2} )^{2} }\) for a bi-objective optimization problem (Since the value of \(Z_{1}\) is much larger than that of \(Z_{2}\), to reflect the influence of \(Z_{2}\) on the MID, \(Z_{1}\) is reduced by 5000 times). Here, the ideal point \((o_{1} ,o_{2} )\) is (0, 0). The lower the value of MID, the better the solution quality.

4.2 Comparative Analysis

For comparing the performance of NSGA-II, EBEGA and MPGA, all algorithms solve 9 test experiments for 10 times and the above performance measures are comparatively analyzed with the same parameters. The experiments are performed on the testing examples in Table 1. By comprehensive tests of the parameters for three evaluation methods, the adopted parameters in experiments 1–3 are: the maximum iteration number \(G = 1000\), the number of chromosomes \({\text{popsize}} = 30\) (\({\text{popsize}} = 15\) of each sub-population for MPGA). For experiments 4–6, \(G = 1200\), \({\text{popsize}} = 40\) (\({\text{popsize}} = 20\) of each sub-population for MPGA). For experiments 7–9, \(G = 1500\), \({\text{popsize}} = 50\)(\({\text{popsize}} = 25\) of each sub-population for MPGA). In addition, the probability of partheno crossover operator \(P_{c1} { = }0.2\), the probability of arithmetic crossover \(P_{c2} { = }0.1\), the probability of production mutation \(P_{m1} { = }0.4\) and the probability of mutation for the number of workers \(P_{m2} { = }0.5\) when the iteration number is less than 600, otherwise, \(P_{c1} { = }0.3\),\(P_{c2} { = }0.2\), \(P_{m1} { = }0.6\) and \(P_{m2} { = }0.7\). The parameters of labors are shown in Table 2.

Table 1 The experiments definition
Table 2 Parameters of labors

The experiments are performed on DELL Vostro 3800-R1846 (Intel Core i3 4130(3.4 GHz 8 GB memory)). The simulated results are reported in Tables 3, 4, 5 and Fig. 9, respectively. The followings show the efficiency of three solution algorithms under each experiment.

  • \({\text{avg}}(Z_{1} )\): EBEGA > MPGA > NSGA-II,

  • \({\text{avg}}(Z_{2} )\): EBEGA > MPGA > NSGA-II,

  • MID: EBEGA > MPGA > NSGA-II,

  • Runtime: EBEGA > NSGA-II > MPGA.

Table 3 Average test results of EBEGA
Table 4 Average test results of MPGA
Table 5 Average test results of NSGA-II
Fig. 9
figure 9

Comparative analysis of performance measures

Both the \({\text{avg}}(Z_{1} )\) and \({\text{avg}}(Z_{2} )\), and also the MID indicate that NSGA-II performs the best. Moreover, the average runtimes of EBEGA and NSGA-II are much smaller than that of MPGA. There is no significant difference in the number of non-dominated sets obtained by the three algorithms.

The \(M_{2}\) in Tables 3, 4, 5 shows that the overall performance of MPGA is much better than that of EBEGA, the performance of NSGA-II is much better than that of MPGA for all of the testing examples.

The obtained solutions of three algorithms are compared in Fig. 10. The results show that the obtained solution set by NSGA-II is superior to that obtained by EBEGA and MPGA. Then, the evaluation method of NSGA-II is used to analyze the effect of local search algorithm, and the results are shown in Fig. 11. It can be seen from the figure that the effect of local search algorithm is better than that of no local search, and the quality of solutions has been significantly improved.

Fig. 10
figure 10figure 10

Comparison of calculation results for three strategies

Fig. 11
figure 11

Effect analysis of local search algorithm

Comparing the effect of experiment 1, 5 and 9, it can be seen that the local search algorithm can improve the quality of the solution obviously and has a better effect in large-scale experiment. Nevertheless, for each of the experiments, the runtime is greatly increased. For example, the runtime will be increased from 15.55 to 93.64 for experiment 1, 30.79–172.40 for experiment 5 and 52.22–319.11 for experiment 9, respectively.

5 Conclusions

The diversity of products and fierce competition make the stability of manufacturing industry more important. To enhance the adaptability of enterprises and reduce production costs, enterprises need to arrange production plan reasonably. Considering the relationship of raw materials, inventory cost and product demand, this paper develops a bi-objective optimization model for an APP problem of multi-product, multi-stage to minimize total production costs and instability in the work force. Based on the analysis of the feasible range of the planned production and the number of workers in each period, a genetic algorithm is designed to solve the model combined with local search algorithm.

Moreover, 9 test experiments are employed to validate the performance of the proposed evaluation of NSGA-II with the ones of EBEGA and MPGA. The performances of these three algorithms were compared statistically. The computational results indicated that the proposed evaluation of NSGA-II with local search algorithm had high computational performance, and the designed local search algorithm could improve the quality of solutions significantly, but it needed extra computing time.

In this paper, we assumed that the product demand and production time are known. In reality, product demand and production time may be uncertain. Therefore, future research is needed to consider the uncertainty to design the related APP model and algorithm.