A Production-Inventory Model for a Deteriorating Item Incorporating Learning Effect Using Genetic Algorithm

Abstract

Demand for a seasonal product persists for a fixed period of time. Normally the “finite time horizon inventory control problems” are formulated for this type of demands. In reality, it is difficult to predict the end of a season precisely. It is thus represented as an uncertain variable and known as random planning horizon. In this paper, we present a production-inventory model for deteriorating items in an imprecise environment characterised by inflation and timed value of money and considering a constant demand. It is assumed that the time horizon of the business period is random in nature and follows exponential distribution with a known mean. Here, we considered the resultant effect of inflation and time value of money as both crisp and fuzzy. For crisp inflation effect, the total expected profit from the planning horizon is maximized using genetic algorithm (GA) to derive optimal decisions. This GA is developed using Roulette wheel selection, arithmetic crossover, and random mutation. On the other hand when the inflation effect is fuzzy, we can expect the profit to be fuzzy, too! As for the fuzzy objective, the optimistic or pessimistic return of the expected total profit is obtained using, respectively, a necessity or possibility measure of the fuzzy event. The GA we have developed uses fuzzy simulation to maximize the optimistic/pessimistic return in getting an optimal decision. We have provided some numerical examples and some sensitivity analyses to illustrate the model.

1. Introduction

Existing theories of inventory control implicitly assumed that lifetime of the product is infinite and models are developed under finite or infinite planning horizon such as that of Bartmann and Beckmann [1], Hadley and Whitin [2], Roy et al. [3], and Roy et al. [4]. In reality, however, products rarely have an infinite lifetime, and there are several reasons for this. Change in product specifications and design may lead to a newer version of the product. Sometimes, due to rapid development of technology (cf. Gurnani [5]), a product may be abandoned, or even be substituted by another product. On the other hand, assuming a finite planning horizon is not appropriate, for example, for a seasonal product, though planning horizon is normally assumed as finite and crisp, it fluctuates in every year depending upon the rate of production, environmental effects, and so forth. Hence, it is better to estimate this horizon as having a fuzzy or stochastic nature. Moon and Yun [6] developed an Economic Ordered Quantity (EOQ) model in a random planning horizon. Moon and Lee [7] further developed an EOQ model taking account of inflation and time discounting, with random product life cycles. Recently, Roy et al. [8] and Roy et al. [9], developed inventory models with stock-dependent demand over a random planning horizon under imprecise inflation and finite discounting. Yet, till now, none has developed an Economic Production Quantity (EPQ) model, which incorporates the lifetime of a product as a random variable.

Production cost of a manufacturing system depends upon the combination of different production factors. These factors are (a) raw materials, (b) technical knowledge, (c) production procedure, (d) firm size, (e) quality of product and so forth, Normally, the cost of raw materials is imprecise in nature. So far, cost of technical knowledge, that is, labor cost, has been usually assumed to be constant. However, because the firms and employees perform the same task repeatedly, they learn how to repeatedly provide a standard level of performance. Therefore, processing cost per unit product decreases in every cycle. Similarly part of the ordering cost may also decrease in every cycle. In the inventory control literature, this phenomenon is known as the learning effect. Although different types of learning effects in various areas have been studied (cf. Chiu and Chen [10], Kuo and Yang [11], Alamri and Balkhi [12], etc.), it has rarely been studied in the context of inventory control problems.

Several studies have examined the effect of inflation on inventory policy. Buzacott [13] first developed an approach on modelling inflation-assuming constant inflation rate subject to different types of pricing policies. Misra [14] proposed an inflation model for the EOQ, in which the time value of money and different inflation rates were considered. Brahmbhatt [15] also developed an EOQ model under a variable inflation rate and marked-up prices. Later, Gupta and Vrat [16] developed a multi-item inventory model for a variable inflation rate. Though a considerable number of researches (cf. Padmanabhan and Vrat [17], Hariga and Ben-Daya [18], Chen [19], Dey et al. [20], etc.) have been done in this area, none has considered the imprecise inflationary effect on EPQ model, especially when the lifetime of the product is random.

In dealing with these shortcomings above, this paper shows an EPQ model of a deteriorating item with a random planning horizon, that is, the lifetime of the product is assumed as random in nature and it follows an exponential distribution with a known mean. Unit production cost decreases in each production cycle due to learning effects of the workers on production. Similarly, setup cost in each cycle is partly constant and partly decreasing in each cycle due to learning effects of the employees. The model is formulated to maximize the expected profit from the whole planning horizon and is solved using genetic algorithm (GA). It is illustrated with some numerical data, and some sensitivity analyses on expected profit function are so presented.

2. Assumptions and Notations

In this paper, the mathematical model is developed on the basis of the following assumptions and notations.

Assumptions
(1)Demand rate is known and constant. (2)Time horizon (a random variable) is finite. (3)Time horizon accommodates first 𝑁 cycles and ends during (𝑁+1) cycles. (4)Setup time is negligible. (5)Production rate is known and constant. (6)Shortages are not allowed. (7)A constant fraction of on-hand inventory gets deteriorated per unit time. (8)Lead time is zero. (9)Production cost and setup cost decrease due to the learning in setups and improvement in quality.

Notations
The notations used in this paper are listed below.𝑞(𝑡): on hand inventory of a cycle at time 𝑡, (𝑗−1)𝑇≤𝑡≤𝑗𝑇 (𝑗=1,2,…,𝑁).𝑡₁: production period in each cycle.𝑃: Production rate in each cycle.𝐷: demand rate in each cycle.𝐶₁: holding cost per unit item per unit time. 𝐶^𝑗₃=𝐶₃+𝐶^′₃𝑒^−𝛽𝑗: is setup cost in 𝑗th(𝑗=1,2,…,𝑁) cycle, 𝛽>0 (𝛽 is the learning coefficient associated with setup cost).𝑝₀𝑒^−𝛾𝑗: production cost in the 𝑗th(𝑗=1,2,…,𝑁) cycle, 𝑝₀,𝛾>0 (𝛾 is the learning coefficient associated with production cost).𝑚₀𝑝₀𝑒^−𝛾𝑗: selling price in the 𝑗th(𝑗=1,2,…,𝑁) cycle, 𝑝₀,𝛾>0, 𝑚₀>1.𝑁: number of fully accommodated cycles to be made during the prescribed time horizon.𝑇: duration of a complete cycle.𝑖: inflation rate.𝑟: discount rate.𝑅: 𝑟-𝑖, may be crisp or fuzzy.𝑃(𝑁,𝑇): total profit after completing 𝑁 fully accommodated cycles.𝐻: total time horizon (a random variable) and ℎ is the real time horizon.𝑚₁𝑝₀𝑒^{−𝛾(𝑁+1)}: reduced selling price for the inventory items in the last cycle at the end of time horizon, 𝑝₀,𝛾> 0, 𝑚₁<1.𝜃: deterioration rate of the produced item.𝐸{𝑃(𝑁,𝑇)}: expected total profit from 𝑁 complete cycles.𝐸{𝑇𝑃_𝐿(𝑇)}: expected total profit from the last cycle.𝐸{𝑇𝑃(𝑇)}: expected total profit from the planning horizon.

3. Mathematical Formulation

In this section, we formulate a production-inventory model for deteriorating items under inflation over a random planning horizon incorporating learning effect. Here we assume that there are 𝑁 full cycles during the real time horizon ℎ and the planning horizon ends within the (𝑁+1)th cycle, that is, within the time 𝑡=𝑁𝑇 and 𝑡=(𝑁+1)𝑇. At the beginning of every 𝑗th (𝑗=1,2,…𝑁+1) cycle production starts at 𝑡=(𝑗−1)𝑇 and continues up to 𝑡=(𝑗−1)𝑇+𝑡₁, inventory gradually increases after meeting the demand due to production (cf. Figures 1(a) and 1(b)). Production thus stops at 𝑡=(𝑗−1)𝑇+𝑡₁, and the inventory falls to zero level at the end of the cycle time 𝑡=𝑗𝑇, due to deterioration and consumption. This cycle repeats again and again. For the last cycle some amount may be left after the end of planning horizon. This amount is sold at a reduced price in a lot.

(a)

(b)

Figure 1 

(a). Inventory level when 𝑁𝑇<ℎ<𝑁𝑇+𝑡1. (b) Inventory level when 𝑁𝑇+𝑡1<ℎ<(𝑁+1)𝑇.

Here, it is assumed that the planning horizon 𝐻 is a random variable and follows exponential distribution with probability density function (p.d.f) as 𝑓(ℎ)=𝜆𝑒^−𝜆ℎ,ℎ≥0,0,otherwise. (3.1)

3.1. Formulation for 𝑁 Full Cycles

The differential equations describing the inventory level 𝑞(𝑡) in the interval (𝑗−1)𝑇≤𝑡≤𝑗𝑇(1≤𝑗≤𝑁), 𝑗=1,2,…,𝑁 are given by𝑑𝑞(𝑡)=𝑑𝑡𝑃−𝐷−𝜃𝑞(𝑡),(𝑗−1)𝑇≤𝑡≤(𝑗−1)𝑇+𝑡₁,−𝐷−𝜃𝑞(𝑡),(𝑗−1)𝑇+𝑡₁≤𝑡≤𝑗𝑇, (3.2) where 𝑃>0, 𝐷>0, 𝜃>0, and 0<𝑡₁<𝑇, subject to the conditions that 𝑞(𝑡)=0 at 𝑡=(𝑗−1)𝑇 and 𝑞(𝑡)=0 at 𝑡=𝑗𝑇.

The solutions of the differential equations (3.2) are given by 𝑞(𝑡)=𝑃−𝐷𝜃[1−𝑒^{𝜃{(𝑗−1)𝑇−𝑡}}],(𝑗−1)𝑇≤𝑡≤(𝑗−1)𝑇+𝑡₁,𝐷𝜃[𝑒^{𝜃(𝑗𝑇−𝑡)}−1],(𝑗−1)𝑇+𝑡₁≤𝑡≤𝑗𝑇. (3.3) Now at 𝑡=(𝑗-1)𝑇+𝑡₁, from (3.3), we get𝑃−𝐷𝜃[1−𝑒^−𝜃𝑡1𝐷]=𝜃[𝑒^{𝜃(𝑇−𝑡1)}−1]⇒𝑡₁=1𝜃𝐷ln[1+𝑃(𝑒^𝜃𝑇−1)]. (3.4)

3.2. Total Expected Profit from 𝑁 Full Cycles

From the symmetry of every full cycle, present value of total expected profit from 𝑁 full cycles, 𝐸{𝑃(𝑁,𝑇)}, is given by𝐸{𝑃(𝑁,𝑇)}=ESRN−EPCN−EHCN−ETOCN. (3.5) where ESRN, EPCN, EHCN, and ETOCN are present value of expected total sales revenue, present value of expected total production cost, present value of expected holding cost, and present value of expected total ordering cost, respectively, from 𝑁 full cycles, and their expressions are derived in Appendix A.1 (see (A.13), (A.7), (A.4), (A.10), resp.).

3.3. Formulation for Last Cycle

Duration of the last cycle is [𝑁𝑇,ℎ], where ℎ is the real time horizon corresponding to the random time horizon 𝐻.

Here two different cases may arise depending upon the cycle length.

Case 1. 𝑁𝑇<ℎ≤𝑁𝑇+𝑡₁.

Case 2. 𝑁𝑇+𝑡₁<ℎ≤(𝑁+1)𝑇.

The differential equation describing the inventory level 𝑞(𝑡) in the interval 𝑁𝑇<𝑡≤ℎ are given by𝑑𝑞(𝑡)=𝑑𝑡𝑃−𝐷−𝜃𝑞(𝑡),𝑁𝑇≤𝑡≤𝑁𝑇+𝑡₁,−𝐷−𝜃𝑞(𝑡),𝑁𝑇+𝑡₁≤𝑡≤(𝑁+1)𝑇. (3.6) subject to the conditions that𝑞(𝑁𝑇)=0,𝑞{(𝑁+1)𝑇}=0. (3.7) The solutions of the differential equations in (3.6) are given by𝑞(𝑡)=𝑃−𝐷𝜃[1−𝑒^{𝜃(𝑁𝑇−𝑡)}],𝑁𝑇≤𝑡≤𝑁𝑇+𝑡₁,𝐷𝜃[𝑒^{𝜃{(𝑁+1)𝑇−𝑡}}−1],𝑁𝑇+𝑡₁≤𝑡≤(𝑁+1)𝑇. (3.8)

3.4. Expected Total Profit from Last Cycle

Present value of expected total profit from last cycle is given by𝐸{𝑇𝑃_𝐿(𝑇)}=ESR_𝐿+ERSP_𝐿−EHC_𝐿−EPC_𝐿−EOC_𝐿. (3.9) where ESR_𝐿,ERSP_𝐿,EHC_𝐿,EPC_𝐿,and EOC_𝐿 are present value of expected sales revenue, present value of expected reduced selling price, present value of expected holding cost, present value of expected production cost, present value of expected ordering cost, respectively, from the last cycle, and their expressions are derived in Appendix A.2 (see (A.24), (A.26), (A.20), (A.23), and (A.25), resp.).

3.5. Total Expected Profit from the System

Now, total expected profit from the complete time horizon is given by𝐸{𝑇𝑃(𝑇)}=𝐸{𝑃(𝑁,𝑇)}+𝐸{𝑇𝑃_𝐿(𝑇)}. (3.10)

4. Problem Formulation

4.1. Stochastic Model (Model-1)

When the resultant effect of inflation and discounting (𝑅) is crisp in nature, then our problem is to determine 𝑇 toMax𝐸(𝑇𝑃),subjectto𝑇≥0. (4.1)

4.2. Fuzzy Stochastic Model (Model-2)

In the real world, resultant effect of inflation and time value of money (𝑅) is imprecise, that is, vaguely defined in some situations. So we take 𝑅 as fuzzy number, denoted by 𝑅. Then, due to this assumption, our objective function 𝐸(𝑇𝑃) becomes 𝐸(𝑇𝑃). Since optimization of a fuzzy objective is not well defined, so instead of 𝐸(𝑇𝑃) one can optimize its equivalent optimistic or pessimistic return of the objective as proposed by M. K. Maiti and M. Maiti [21]. Using this method the problem can be reduced to an equivalent crisp problem as discussed below.

If 𝐴 and 𝐵 are two fuzzy subsets of real numbers ℜ with membership functions 𝜇_𝐴 and 𝜇_𝐵, respectively, then taking degree of uncertainty as the semantics of fuzzy number, according to Liu and Iwamura [22], Dubois and Prade [23, 24], and Zimmermann [25],Pos(𝐴⋆𝐵)=sup{min(𝜇_𝐴(𝑥),𝜇_𝐵(𝑦)),𝑥,𝑦∈ℜ,𝑥⋆𝑦}, (4.2) where the abbreviation Pos represent possibility and ⋆ is any one of the relations >,<,=,≤,≥.

On the other hand necessity measure of an event 𝐴⋆𝐵 is a dual of possibility measure. The grade of necessity of an event is the grade of impossibility of the opposite event and is defined asNes(𝐴⋆𝐵)=1−Pos(𝐴⋆𝐵), (4.3) where the abbreviation Nes represents necessity measure and 𝐴⋆𝐵 represents complement of the event 𝐴⋆𝐵.

So for the fuzzy stochastic model one can maximize the crisp variable 𝑧 such that necessity/possibility measure of the event {𝐸(𝑇𝑃)>𝑧} exceeds some predefined level according to decision maker in pessimistic/optimistic sense. Accordingly the problem reduces to the following two models.

Model-2a
When the decision maker prefers to optimize the optimistic equivalent of 𝐸(𝑇𝑃), the problem reduces to determine 𝑇 to Maximize𝑧subjectetopos{𝐸(𝑇𝑃)≥𝑧}≥𝛼₁, (4.4) where 𝛼₁ is confidence level.

Model-2b
On the other hand when the decision maker desires to optimize the pessimistic equivalent of 𝐸(𝑇𝑃), the problem is reduced to determine 𝑇 to Maximize𝑧subjectto,nes{𝐸(𝑇𝑃)≥𝑧}≥𝛼₂thatis,pos{𝐸(𝑇𝑃)≤𝑧}<1−𝛼₂, (4.5) where 𝛼₂ is confidence level.

5. Solution Methodology

To solve the stochastic model (model-1), genetic algorithm (GA) and simulated annealing (SA) are used. The basic technique to deal with problem (4.4) or (4.5) is to convert the possibility/necessity constraint to its deterministic equivalent. However, the procedure is usually very hard and successful in some particular cases (cf. M. K. Maiti and M. Maiti [21]). Following Liu and Iwamura [22] and M. K. Maiti and M. Maiti [21], here two simulation algorithms are proposed to determine 𝑧 in (4.4) and (4.5), respectively, for a feasible 𝑇.

Algorithm 1. Algorithm to determine a feasible 𝑇 to evaluate 𝑧 for the problem(4.4)
To determine 𝑧 for a feasible 𝑇, roughly find a point 𝑅₀ from fuzzy number 𝑅, which approximately minimizes 𝑧. Let this value be 𝑧₀ and set 𝑧=𝑧₀ (For simplicity one can take 𝑧₀=0). Then 𝑅₀ is randomly generated in 𝛼₁-cut set of 𝑅 and let 𝑧₀= value of 𝐸(𝑇𝑃) for 𝑅=𝑅₀ and if 𝑧<𝑧₀ replace 𝑧 with 𝑧₀. This step is repeated a finite number of times and final value is taken as the value of 𝑧. This phenomenon is used to develop the algorithm.(1)Set 𝑧=𝑧₀.(2)Generate 𝑅₀ uniformly from the 𝛼₁ cut set of fuzzy number 𝑅. (3)Set 𝑧₀= value of 𝐸(𝑇𝑃) for 𝑅=𝑅₀.(4)If 𝑧<𝑧₀ then set 𝑧=𝑧₀.(5)Repeat steps 2, 3 and 4, 𝑁₁ times, where 𝑁₁ is a sufficiently large positive integer. (6)Return 𝑧.(7)End algorithm.

Algorithm 2. Algorithm to determine a feasible T to evaluate 𝑧 for the problem (4.5):
We know that nes{𝐸(𝑇𝑃)≥𝑧}≥𝛼₂⇒pos{𝐸(𝑇𝑃)<𝑧}≤1−𝛼₂. Now roughly find a point 𝑅₀ from fuzzy number 𝑅, which approximately minimizes 𝐸(𝑇𝑃). Let this value be 𝑧₀ (For simplicity one can take 𝑧₀=0 also) and 𝜀 be a positive number. Set 𝑧=𝑧₀−𝜀 and if pos{𝐸(𝑇𝑃)<𝑧}≤1−𝛼₂ then increase 𝑧 with 𝜀. Again check pos{𝐸(𝑇𝑃)<𝑧}≤1−𝛼₂ and it continues until 𝑝𝑜𝑠{𝐸(𝑇𝑃)<𝑧}>1−𝛼₂. At this stage decrease value of 𝜀 and again try to improve 𝑧. When 𝜀 becomes sufficiently small then we stop and final value of 𝑧 is taken as the value of 𝑧. Using this criterion, required algorithm is developed as below. In the algorithm the variable 𝐹₀ is used to store initial assumed value of 𝑧 and 𝐹 is used to store value of 𝑧 in each iteration. (1)Set 𝑧=𝑧₀−𝜀, 𝐹=𝑧₀−𝜀, 𝐹₀=𝑧₀−𝜀, tol=0.0001. (2)Generate 𝑅₀ uniformly from the 1−𝛼₂ cut set of fuzzy number 𝑅. (3)Set 𝑧₀= value of 𝐸(𝑇𝑃) for 𝑅=𝑅₀. (4)If 𝑧₀<𝑧.(5)   Then go to step 11.(6)End If(7)Repeat step-2 to step-6 𝑁₂ times.(8)Set 𝐹=𝑧.(9)Set 𝑧=𝑧+𝜀.(10)Go to step-2.(11)If (𝑧=𝐹₀) // In this case optimum value of 𝑧<𝑧₀−𝜀(12)   Set 𝑧=𝐹₀−𝜀, 𝐹=𝐹−𝜀, 𝐹₀=𝐹₀−𝜀. (13)   Go to step-2(14)End If(15)If (𝜀<tol)(16)   Go to step-21(17)End If(18)𝜀=𝜀/10(19)𝑧=𝐹+𝜀(20)Go to step-2.(21)Output 𝐹.(22)End algorithm.

So for a feasible value of 𝑇, we determine 𝑧 using the above algorithms, and to optimize 𝑧 we use GA. GA used to solve model-1 is presented below. When fuzzy simulation algorithm is used to determine 𝑧 in the algorithm, this GA is named fuzzy simulation-based genetic algorithm (FSGA). This is used to determine fuzzy objective function values.

5.1. Genetic Algorithm (GA)/Fuzzy Simulation-Based Genetic Algorithm (FSGA)

Genetic Algorithm is a class of adaptive search technique based on the principle of population genetics. In natural genesis. we know that chromosomes are the main carriers of the hereditary information from parents to offsprings and that genes, which carry hereditary factors, are lined up in chromosomes. At the time of reproduction, crossover and mutation take place among the chromosomes of parents. In this way, hereditary factors of parents are mixed up and carried over to their offsprings. Darwinian principle states that only the fittest animals can survive in nature. So a pair of the fittest parents normally reproduce better offspring.

The above- mentioned phenomenon is followed to create a genetic algorithm for an optimization problem. Here potential solutions of the problem are analogous with the chromosomes and chromosome of better offspring with the better solution of the problem. Crossover and mutation are performed among a set of potential solutions, and a new set of solutions are obtained. It continues until terminating conditions are encountered. Michalewicz [26] proposed a genetic algorithm named the Contractive Mapping Genetic Algorithm (CMGA) and proved the asymptotic convergence of the algorithm by the Banach fixed-point theorem. In CMGA, movement from an old population to a new population takes place only when the average fitness of a new population is better than the old one. This algorithm is modified with the help of a fuzzy simulation process to solve the fuzzy stochastic models of this paper. The algorithm is named FSGA, and this is presented below. In the algorithm, 𝑝_𝑐,𝑝_𝑚 are probabilities of the crossover and the probability of mutation, respectively, 𝐼 is the iteration counter, and 𝑃(𝐼) is the population of potential solutions for iteration 𝐼. The (𝑃(𝐼)) function initializes the population 𝑃(𝐼) at the time of initialization. The (𝑃(𝐼)) function evaluates the fitness of each member of 𝑃(𝐼), and at this stage an objective function value due to each solution is evaluated via the fuzzy simulation process (using algorithm 1 or algorithm 2). In case of stochastic model (model-1) objective function is evaluated directly without using simulation algorithms. So in that case this GA is named ordinary GA. 𝑀 is iteration counter in each generation to improve 𝑃(𝐼), and 𝑀₀ is upper limit of 𝑀.

5.2. GA/FSGA Algorithm

(1)Set 𝐼=0, 𝑀=0, 𝑀₀=50.(2)Initialize 𝑝_𝑐,𝑝_𝑚.(3)Initialize (𝑃(𝐼)) and let 𝑁^′ be its size.(4)Evaluate (𝑃(𝐼)).(5)While (𝑀<𝑀₀)(6) Select 𝑁^′ solutions from 𝑃(𝐼) for mating pool using roulette-wheel selection process [Michalewicz [26]]. Let this set be 𝑃₁(𝐼).(7) Select solutions from 𝑃₁(𝐼) for crossover depending on 𝑝_𝑐.(8) Perform crossover on selected solutions to obtain population 𝑃₁(𝐼).(9) Select solutions from 𝑃₁(𝐼) for mutation depending on 𝑝_𝑚.(10) Perform mutation on selected solutions to obtain new population 𝑃(𝐼+1).(11) Evaluate (𝑃(𝐼+1)).(12)  Set 𝑀=𝑀+1.(13) If average fitness of 𝑃(𝐼+1)> average fitness of 𝑃(𝐼) then(14)  Set 𝐼=𝐼+1.(15)  Set 𝑀=0.(16) End If.(17)End While.(18)Output: Best solution of 𝑃(𝐼).(19)End algorithm.

5.3. GA/FSGA Procedures

(a) Representation
An “𝑛 dimensional real vector” 𝑋=(𝑥₁,𝑥₂,…,𝑥_𝑛) is used to represent a solution, where 𝑥₁, 𝑥₂,…, 𝑥_𝑛 represent 𝑛 decision variables of the problem.

(b) Initialization
𝑁^′ such solutions 𝑋₁, 𝑋₂, 𝑋₃,…, 𝑋_𝑁′ are randomly generated by random number generator. This solution set is taken as initial population 𝑃(𝐼). Here we take 𝑁^′=50, 𝑝_𝑐=0.3, 𝑝_𝑚=0.2, and 𝐼=1. These parametric values are assumed as these giving better convergence of the algorithm for the model.

(c) Fitness value
Value of the objective function due to the solution 𝑋 is taken as fitness of 𝑋. Let it be 𝑓(𝑋). Objective function is evaluated via fuzzy simulation process (using Algorithm 1 or Algorithm 2) for model-2.

(d) Selection Process for Mating Pool
The following steps are followed for this purpose.(i)Find total fitness of the population 𝐹=∑^𝑁′_𝑖=1𝑓(𝑋_𝑖). (ii)Calculate the probability of selection 𝑝_𝑖 of each solution 𝑋_𝑖 by the formula 𝑝_𝑖=𝑓(𝑋_𝑖)/𝐹. (iii)Calculate the cumulative probability 𝑞_𝑖 for each solution 𝑋_𝑖 by the formula 𝑞_𝑖=∑^𝑖_𝑗=1𝑝_𝑗. (iv)Generate a random number “𝑟” from the range [0,1]. (v)If 𝑟<𝑞₁, then select 𝑋₁: otherwise select 𝑋_𝑖(2≤𝑖≤𝑁), where 𝑞_𝑖−1≤𝑟≤𝑞_𝑖.(vi) Repeat step (iv) and (v) 𝑁^′ times to select 𝑁^′ solutions from old population. Clearly one solution may be selected more than once. (vii)Selected solution set is denoted by 𝑃₁(𝐼) in the proposed GA/FSGA algorithm.

(c) Crossover
(i)Selection for Crossover. For each solution of 𝑃(𝐼) generate a random number 𝑟 from the range [0,1]. If 𝑟<𝑝_𝑐, then the solution is taken for crossover, where 𝑝_𝑐 is the probability of crossover. (ii)Crossover Process. Crossover takes place on the selected solutions. For each pair of coupled solutions 𝑌₁, 𝑌₂, a random number 𝑐 is generated from the range [0,1] and their offsprings 𝑌₁₁ and 𝑌₂₁ are obtained by the formula 𝑌₁₁=𝑐𝑌₁+(1−𝑐)𝑌₂,𝑌₂₁=𝑐𝑌₂+(1−𝑐)𝑌₁. (5.1)

(d) Mutation
(i)Selection for Mutation. For each solution of 𝑃(𝐼) generate a random number 𝑟 from the range [0,1]. If 𝑟<𝑝_𝑚, then the solution is taken for mutation, where 𝑝_𝑚 is the probability of mutation. (ii)Mutation Process. To mutate a solution 𝑋=(𝑥₁,𝑥₂,…,𝑥_𝑛) select a random integer 𝑟 in the range [1, 𝑛]. Then replace 𝑥_𝑟 by randomly generated value within the boundary of the 𝑟th component of 𝑋.

6. Numerical Illustration

6.1. Stochastic Model

The following numerical data are used to illustrate the model:

𝐶₃=$50,𝐶^′₃=$100,𝐶₁=$0.75,𝛾=0.05,𝛽=0.5, 𝜆=0.01,𝑚₀=1.8,𝑚₁=0.8,𝑟=0.1,𝑖=0.05, that is 𝑅=0.05, 𝜃=0.1, 𝑝₀=4 in appropriate units.

The fuzzy simulation-based GA designed in Section 5.3 is used to solve the model. Here, the initial population size is 50, the probability of crossover is 0.3, and the probability of mutation is 0.2. After 50 iterations the results obtain are shown in Table 1(a). The optimal values of 𝑇 along with maximum expected total profit have been calculated for different values of 𝑃 and 𝐷, and results in GA are displayed in Table 1(a). In order to verify the feasibility of our proposed algorithm we combine a Simulated Annealing (Appendix B) to solve the same numerical example. The result using SA is displayed in Table 1(b).

Table 1 

Comparison of Results Using GA and SA
It is observed that in all cases genetic algorithm (GA) gives the better results than simulated annealing (SA). Also it is observed that in GA after fifty iterations we get the above results but in SA we get the results by taking more than fifty iterations. Accordingly, the performance of GA is acceptable.

Sensitivity Analysis
Sensitivity analysis is performed for stochastic model with respect to different 𝜆, 𝛽, 𝛾, and 𝑅 values for crisp inflation, and results are presented in Tables 2, 3, 4, and 5, and Figures 2, 3, 4, and 5, respectively, when other input values are the same. It is observed that profit decreases and 𝜆 increases; when 𝛽 increases, setup cost decreases and as such profit increases; also when 𝛾 increases, unit production cost(𝑝₀) decreases, as well as selling price also decreases, then profit decreases and profit decreases with R increases, which agrees with reality.

Table 2 

Results due to different 𝜆.

Table 3 

Results due to different 𝛽.

Table 4 

Results due to different 𝛾.

Table 5 

Results due to different 𝑅.

Figure 2 

Figure 3 

Figure 4 

Figure 5 

6.2. Fuzzy Stochastic Model

Here the resultant inflationary effect is considered as a triangular fuzzy number, that is, 𝑅 = 𝑟̃−𝑖̃ = (0.095, 0.1, 0.105) − (0.045, 0.05, 0.055) = (0.04, 0.05, 0.06), and all other data remain the same as in stochastic model. The maximum optimistic/pessimistic return from expression (4.4), (4.5) has been calculated for different values of possibility and necessity, and results are displayed in Table 6.

Table 6 

Results due to possibility and necessity.

7. Conclusion

In this paper, for the first time an economic production quantity model for deteriorating items has been considered under inflation and time discounting over a stochastic time horizon. Also for the first time learning effect on production and setup cost is incorporated in an economic production quantity model. The methodology presented here is quite general and provides a valuable reference for decision makers in the production inventory system. To solve the proposed highly nonlinear models, we have designed a fuzzy simulation based GA. The algorithm has been tested using a numerical example. The results show that the algorithms designed in the paper perform well. Finally, a future study will incorporate more realistic assumptions in the proposed model, such as variable demand and production, allowing shortages and so forth.

Appendices

A.

A.1. Calculation for Expected Sales Revenue for 𝑁 Full Cycles

Present value of holding cost of the inventory for the 𝑗th(1≤𝑗≤𝑁) cycle, (𝐻𝐶_𝑗), is given by𝐻𝐶_𝑗=𝐶₁^{(𝑗−1)𝑇+𝑡1}_{(𝑗−1)𝑇}𝑞(𝑡)𝑒^−𝑅𝑡𝑑𝑡+𝐶₁^𝑗𝑇_{(𝑗−1)𝑇+𝑡1}𝑞(𝑡)𝑒^−𝑅𝑡=𝐶𝑑𝑡₁(𝑃−𝐷)𝜃𝑅[𝑒^{−𝑅(𝑗−1)𝑇}−𝑒^{−𝑅{(𝑗−1)𝑇+𝑡1}}]−𝐶₁(𝑃−𝐷)𝜃(𝜃+𝑅)[𝑒^{−𝑅(𝑗−1)𝑇}−𝑒^{−𝑅{(𝑗−1)𝑇+𝑡1}−𝜃𝑡1}]+𝐶₁𝐷𝜃(𝜃+𝑅)[𝑒^{𝜃𝑗𝑇−(𝜃+𝑅){(𝑗−1)𝑇+𝑡1}}−𝑒^{−𝑅𝑗𝑇}]+𝐶₁𝐷𝜃𝑅[𝑒^{−𝑅𝑗𝑇}−𝑒^{−𝑅{(𝑗−1)𝑇+𝑡1}}].Also,^𝑁_𝑗=1𝑒^{−𝑅(𝑗−1)𝑇}=(1−𝑒^{−𝑁𝑅𝑇}1−𝑒^−𝑅𝑇). (A.1)(A.2)

Total holding cost from 𝑁 full cycles, (HCN), is given byHCN=^𝑁_𝑗=1𝐻𝐶_𝑗𝐶=[¹(𝑃−𝐷)𝜃𝑅(1−𝑒^−𝑅𝑡1𝐶)−¹(𝑃−𝐷)𝜃(𝜃+𝑅)(1−𝑒^{−(𝑅+𝜃)𝑡1})−𝐶₁𝐷𝜃(𝜃+𝑅)(1−𝑒^{(𝜃+𝑅)(𝑇−𝑡1)})𝑒^−𝑅𝑇+𝐶₁𝐷𝜃𝑅(1−𝑒^{𝑅(𝑇−𝑡1)})𝑒^−𝑅𝑇](1−𝑒^{−𝑁𝑅𝑇}1−𝑒^−𝑅𝑇). (A.3) So, the present value of expected holding cost from 𝑁 complete cycles, (EHCN), is given byEHCN=^∞_𝑁=0^(𝑁+1)𝑇_𝑁𝑇𝐶HCN.𝑓(ℎ)𝑑ℎ=[¹(𝑃−𝐷)(𝜃𝑅1−𝑒^−𝑅𝑡11−𝑒^−𝑅𝑇𝐶)−¹(𝑃−𝐷)(𝜃(𝜃+𝑅)1−𝑒^{−(𝑅+𝜃)𝑡1}1−𝑒^−𝑅𝑇)−𝐶₁𝐷(𝜃(𝜃+𝑅)1−𝑒^{(𝜃+𝑅)(𝑇−𝑡1)}1−𝑒^−𝑅𝑇)𝑒^−𝑅𝑇+𝐶₁𝐷(𝜃𝑅1−𝑒^{𝑅(𝑇−𝑡1)}1−𝑒^−𝑅𝑇)𝑒^−𝑅𝑇](1−1−𝑒^−𝜆𝑇1−𝑒^{−(𝑅+𝜆)𝑇}). (A.4) Present value of production cost for the 𝑗th(1≤𝑗≤𝑁) cycle, (𝑃𝐶_𝑗), is given by𝑃𝐶_𝑗=𝑝₀𝑒^−𝛾𝑗⋅𝑃^{(𝑗−1)𝑇+𝑡1}_{(𝑗−1)𝑇}𝑒^−𝑅𝑡𝑝𝑑𝑡=⁰𝑒^−𝛾𝑗⋅𝑃𝑅(1−𝑒^−𝑅𝑡1)𝑒^{−𝑅(𝑗−1)𝑇}. (A.5) Present value of total production cost from 𝑁 full cycles, (PCN), is given byPCN=^𝑁_𝑗=1𝑃𝐶_𝑗=𝑝₀𝑅⋅𝑃⋅𝑒^𝑅𝑇⋅(1−𝑒^−𝑅t1)⋅𝑒^{−(𝛾+𝑅𝑇)}⋅(1−𝑒^{−𝑁(𝛾+𝑅𝑇)}1−𝑒^{−(𝛾+𝑅𝑇)}). (A.6) Present value of expected total production cost from 𝑁 full cycles, (EPCN), is given byEPCN=^∞_𝑁=0^(𝑁+1)𝑇_𝑁𝑇=𝑝PCN⋅𝑓(ℎ)𝑑ℎ₀𝑅⋅𝑃⋅𝑒^𝑅𝑇⋅(1−𝑒^−𝑅𝑡1)⋅𝑒^{−(𝛾+𝑅𝑇)}𝑒⋅{^−𝜆𝑇(1−𝑒^{−(𝛾+𝑅𝑇+𝜆𝑇)})}. (A.7) Present value of ordering cost for the 𝑗th(1≤𝑗≤𝑁) cycle, 𝐶^𝑗₃, is given by𝐶^𝑗₃={𝐶₃+𝐶^′₃⋅𝑒^−𝛽𝑗}⋅𝑒^{−𝑅(𝑗−1)𝑇},𝐶₃,𝐶^′₃,𝛽>0. (A.8) Present value of total ordering cost from 𝑁 full cycles, (TOCN), is given byTOCN=^𝑁_𝑗=1𝐶^𝑗₃=𝐶₃(1−𝑒^{−𝑁𝑅𝑇}1−𝑒^−𝑅𝑇)+𝐶^′₃⋅𝑒^−𝛽⋅(1−𝑒^{−𝑁(𝛽+𝑅𝑇)}1−𝑒^{−(𝛽+𝑅𝑇)}). (A.9) Present value of expected total ordering cost from 𝑁 full cycles, (ETOCN), is given byETOCN=^∞_𝑁=0^(𝑁+1)𝑇_𝑁𝑇=𝐶TOCN⋅𝑓(ℎ)𝑑ℎ₃𝑒^−𝜆𝑇(1−𝑒^{−(𝜆+𝑅)𝑇})+𝐶^′₃⋅𝑒^−𝛽⋅𝑒^−𝜆𝑇(1−𝑒^{−(𝛽+𝑅𝑇+𝜆𝑇)}). (A.10) Present value of sales revenue for the 𝑗th(1≤𝑗≤𝑁) cycle, (𝑆𝑅_𝑗), is given by𝑆𝑅_𝑗=𝑚₀⋅𝑝₀⋅𝑒^{−𝛾𝑗𝑗𝑇}_{(𝑗−1)𝑇}𝐷⋅𝑒^−𝑅𝑡=𝑚𝑑𝑡₀⋅𝑝₀⋅𝑒^−𝛾𝑗⋅𝐷𝑅{𝑒^{−𝑅(𝑗−1)𝑇}−𝑒^{−𝑅𝑗𝑇}}. (A.11) Present value of total sales revenue from 𝑁 full cycles, (SRN), is given bySRN=^𝑁_𝑗=1𝑆𝑅_𝑗=𝑚₀⋅𝑝₀𝑅⋅𝐷⋅(𝑒^𝑅𝑇−1)⋅𝑒^{−(𝛾+𝑅𝑇)}⋅(1−𝑒^{−𝑁(𝛾+𝑅𝑇)}1−𝑒^{−(𝛾+𝑅𝑇)}). (A.12) Present value of expected total sales revenue from 𝑁 full cycles, (ESRN), is given byESRN=^∞_𝑁=0^(𝑁+1)𝑇_𝑁𝑇SRN⋅𝑓(ℎ)𝑑ℎ=𝑚₀⋅𝑝₀𝑅⋅𝐷⋅(𝑒^𝑅𝑇−1)⋅𝑒^{−(𝛾+𝑅𝑇)}𝑒⋅{^−𝜆𝑇(1−𝑒^{−(𝛾+𝑅𝑇+𝜆𝑇)})}. (A.13)

A.2. Calculation for Expected Sales Revenue for Last Cycle

Case 1 (𝑁𝑇<ℎ≤𝑁𝑇+𝑡₁). Present value of holding cost of the inventory for the last cycle is given by 𝐻𝐶_𝐿1=𝐶₁^ℎ_𝑁𝑇𝑞(𝑡)𝑒^−𝑅𝑡=𝐶𝑑𝑡₁(𝑃−𝐷)𝜃[1𝑅(𝑒^{−𝑁𝑅𝑇}−𝑒^−𝑅ℎ1)+𝜃+𝑅{𝑒^{(𝜃𝑁𝑇−𝜃ℎ−𝑅ℎ)}−𝑒^{−𝑁𝑅𝑇}}]. (A.14) Present value of production cost is given by 𝑃𝐶_𝐿1=𝑝₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝑃^ℎ_𝑁𝑇𝑒^−𝑅𝑡=𝑝𝑑𝑡₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝑃𝑅[𝑒^{−𝑅𝑁𝑇}−𝑒^−𝑅ℎ]. (A.15) Present value of ordering cost is given by {𝐶₃+𝐶^′₃⋅𝑒^{−𝛽(𝑁+1)}}𝑒^{−𝑁𝑅𝑇}.
Present value of sales revenue is given by 𝑆𝑅_𝐿1=𝑚₀⋅𝑝₀⋅e^{−𝛾(𝑁+1)}⋅𝐷^ℎ_𝑁𝑇𝑒^−𝑅𝑡=𝑚𝑑𝑡₀⋅𝑝₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝐷𝑅[𝑒^{−𝑅𝑁𝑇}−𝑒^−𝑅ℎ]. (A.16)

Case 2 (𝑁𝑇+𝑡₁<ℎ≤(𝑁+1)𝑇). Present value of holding cost of the inventory for the last cycle is given by 𝐻𝐶_𝐿2=𝐶₁^{𝑁𝑇+𝑡1}_𝑁𝑇𝑞(𝑡)𝑒^−𝑅𝑡𝑑𝑡+𝐶₁^ℎ_{𝑁𝑇+𝑡1}𝑞(𝑡)𝑒^−𝑅𝑡=𝐶𝑑𝑡₁(𝑃−𝐷)𝜃[1𝑅{𝑒^{−𝑁𝑅𝑇}−𝑒^{−𝑅(𝑁𝑇+𝑡1)}𝑒}+^𝜃𝑁𝑇𝜃+𝑅{𝑒^{−(𝜃+𝑅)(𝑁𝑇+𝑡1)}−𝑒^{−(𝜃+𝑅)𝑁𝑇}𝐶}]+¹𝐷𝜃1×[𝑒𝜃+𝑅^{𝜃(𝑁+1)𝑇}{𝑒^{−(𝜃+𝑅)(𝑁𝑇+𝑡1)}−𝑒^{−(𝜃+𝑅)ℎ}1}+𝑅{𝑒^−𝑅ℎ−𝑒^{−𝑅(𝑁𝑇+𝑡1)}}]. (A.17) Present value of production cost is given by 𝑃𝐶_𝐿2=𝑝₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝑃^{𝑁𝑇+𝑡1}_𝑁𝑇𝑒^−𝑅𝑡=𝑝𝑑𝑡₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝑃𝑅[𝑒^{−𝑅𝑁𝑇}−𝑒^{−𝑅(𝑁𝑇+𝑡1)}]. (A.18) Present value of ordering cost is given by {𝐶₃+𝐶^′₃⋅𝑒^{−𝛽(𝑁+1)}}⋅𝑒^{−𝑁𝑅𝑇}.
Present value of sales revenue is given by 𝑆𝑅_𝐿2=𝑚₀⋅𝑝₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝐷^{𝑁𝑇+𝑡1}_𝑁𝑇𝑒^−𝑅𝑡𝑑𝑡+𝑚₀⋅𝑝₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝐷^ℎ_{𝑁𝑇+𝑡1}𝑒^−𝑅𝑡=𝑚𝑑𝑡₀⋅𝑝₀⋅𝑒^{−𝛾(𝑁+1)}⋅𝐷𝑅[𝑒^{−𝑅𝑁𝑇}−𝑒^−𝑅ℎ]. (A.19) Present value of expected holding cost for the last cycle is given by EHC_𝐿=^∞_𝑁=0^(𝑁+1)𝑇_𝑁𝑇𝐻𝐶_𝐿=⋅𝑓(ℎ)𝑑ℎ_{∞𝑁=0𝑁𝑇+𝑡1𝑁𝑇}𝐻𝐶_𝐿1⋅𝑓(ℎ)𝑑ℎ+^∞_𝑁=0^(𝑁+1)𝑇_{𝑁𝑇+𝑡1}𝐻𝐶_𝐿2⋅𝑓(ℎ)𝑑ℎ=EHC_𝐿1+EHC_𝐿2, (A.20) where EHC_𝐿1=𝐶₁(𝑃−𝐷)𝜃[1𝑅(1−𝑒^−𝜆𝑡1𝜆)−𝑅(𝑅+𝜆){1−𝑒^{−(𝜆+𝑅)𝑡1}}+𝜆(𝜃+𝑅)(𝜃+R+𝜆){1−𝑒^{−(𝜃+𝑅+𝜆)𝑡1}}−1𝜃+𝑅(1−𝑒^−𝜆𝑡11)]1−𝑒^{−(𝜆+𝑅)𝑇},EHC_𝐿2=𝐶₁(𝑃−𝐷)𝜃[1𝑅(𝑒^−𝑅𝑡1−1)(𝑒^−𝜆𝑇−𝑒^−𝜆𝑡1)+1𝜃+𝑅{1−𝑒^{−(𝜃+𝑅)𝑡1}}(𝑒^−𝜆𝑇−𝑒^−𝜆𝑡11)]1−𝑒^{−(𝜆+𝑅)𝑇}+𝐶₁𝐷𝜃[1𝜃+𝑅(𝑒^−𝜆𝑡1−𝑒^−𝜆𝑇)𝑒^{𝜃𝑇−(𝜃+𝑅)𝑡1}+𝜆𝑒^𝜃𝑇(𝜃+𝑅)(𝜃+𝑅+𝜆){𝑒^{−(𝜃+𝑅+𝜆)𝑇}−𝑒^{−(𝜃+𝑅+𝜆)𝑡1}}−𝜆𝑅(𝑅+𝜆){𝑒^{−(𝑅+𝜆)𝑇}−𝑒^{−(𝑅+𝜆)𝑡1}}+1𝑅(𝑒^−𝜆𝑇−𝑒^−𝜆𝑡1)𝑒^−𝑅𝑡1]11−𝑒^{−(𝜆+𝑅)𝑇}. (A.21)(A.22) Present value of expected production cost for the last cycle is given by EPC_𝐿=^∞_𝑁=0^(𝑁+1)T_𝑁𝑇𝑃𝐶_𝐿=⋅𝑓(ℎ)𝑑ℎ_{∞𝑁=0𝑁𝑇+𝑡1𝑁𝑇}𝑃𝐶_𝐿1⋅𝑓(ℎ)𝑑ℎ+^∞_𝑁=0^(𝑁+1)𝑇_{𝑁𝑇+𝑡1}𝑃𝐶_𝐿2=𝑝⋅𝑓(ℎ)𝑑ℎ₀⋅𝑒^−𝛾⋅𝑃𝑅[(1−𝑒^−𝜆𝑡11)⋅(1−𝑒^{−(𝑅𝑇+𝜆𝑇+𝛾)})+𝜆{𝑒(𝑅+𝜆)^{−(𝑅+𝜆)𝑡1}−1(1−𝑒^{−(𝑅𝑇+𝜆𝑇+𝛾)})+𝑝}]₀⋅𝑒^−𝛾⋅𝑃𝑅[(1−𝑒^−𝑅𝑡1)⋅(𝑒^−𝜆𝑡1−𝑒^−𝜆𝑇1)⋅(1−𝑒^{−(𝑅𝑇+𝜆𝑇+𝛾)})]. (A.23) Present value of expected sales revenue from the last cycle is given by ESR_𝐿=^∞_𝑁=0^(𝑁+1)𝑇_𝑁𝑇𝑆𝑅_𝐿=⋅𝑓(ℎ)𝑑ℎ_{∞𝑁=0𝑁𝑇+𝑡1𝑁𝑇}𝑆𝑅_𝐿1⋅𝑓(ℎ)𝑑ℎ+^∞_𝑁=0^(𝑁+1)𝑇_{𝑁𝑇+𝑡1}𝑆𝑅_𝐿2=𝑚⋅𝑓(ℎ)𝑑ℎ₀⋅𝑝₀⋅𝐷𝑅⋅𝑒^−𝛾⋅[(1−𝑒^−𝜆𝑇1)⋅(1−𝑒^{−(𝑅𝑇+𝜆𝑇+𝛾)})+𝜆(𝑅+𝜆)(𝑒^{−(𝑅+𝜆)𝑇}1−1)⋅(1−𝑒^{−(𝑅𝑇+𝜆𝑇+𝛾)})]. (A.24) Present value of expected ordering cost for the last cycle is given by EOC_𝐿=^∞_𝑁=0^(𝑁+1)𝑇_𝑁𝑇{𝐶₃+𝐶^′₃⋅𝑒^{−𝛽(𝑁+1)}}⋅𝑒^{−𝑁𝑅𝑇}𝑓(ℎ)𝑑ℎ=𝐶₃(1−𝑒^−𝜆𝑇)(1−𝑒^{−(𝜆+𝑅)𝑇})+𝐶^′₃⋅𝑒^−𝛽⋅(1−𝑒^−𝜆𝑇)(1−𝑒^{−(𝛽+𝜆𝑇+𝑅𝑇)}). (A.25) Present value of expected reduced selling price from the last cycle is given by ERSP_𝐿=𝑚₁𝑝₀^∞_𝑁=0𝑒^{−𝛾(𝑁+1)(𝑁+1)𝑇}_𝑁𝑇𝑒^−𝑅ℎ𝑞(ℎ)⋅𝑓(ℎ)𝑑ℎ=𝑚₁𝑝₀𝑒^∞−𝛾_𝑁=0𝑒^{−𝛾𝑁𝑁𝑇+𝑡1}_𝑁𝑇𝑒^−𝑅ℎ𝑞(ℎ)⋅𝑓(ℎ)𝑑ℎ+𝑚₁𝑝₀𝑒^∞−𝛾_𝑁=0𝑒^{−𝛾𝑁(𝑁+1)𝑇}_{𝑁𝑇+𝑡1}𝑒^−𝑅ℎ𝑞(ℎ)⋅𝑓(ℎ)𝑑ℎ=ERSP_𝐿1+ERSP_𝐿2, (A.26) where ERSP_𝐿1=𝑚₁𝑝₀𝑒^−𝛾𝜆(𝑃−𝐷)𝜃[1𝑅+𝜆{1−𝑒^{−(𝑅+𝜆)𝑡1}}−1𝑅+𝜃+𝜆{1−𝑒^{−(𝑅+𝜆+𝜃)𝑡1}1}]1−𝑒^{−(𝛾+𝑅𝑇+𝜆𝑇)},ERSP_𝐿2=𝑚₁𝑝₀𝑒^−𝛾𝜆𝐷𝜃[1𝑅+𝜆+𝜃{𝑒^{−(𝑅+𝜆+𝜃)𝑡1}−𝑒^{−(𝑅+𝜆+𝜃)𝑇}}𝑒^𝜃𝑇+1𝑅+𝜆{𝑒^{−(𝑅+𝜆)𝑇}−𝑒^{−(𝑅+𝜆)𝑡1}1}]1−𝑒^{−(𝛾+𝑅𝑇+𝜆𝑇)}. (A.27)(A.28)

B. Simulated Annealing

SA is a stochastic search algorithm developed by mimicking the physical process of evolution of a solid in a heat bath to thermal equilibrium. In the early 1980s Kirkpatrick et al. [27, 28] and independently Cerny [29] introduced the concept of annealing in optimization.

Consider an ensemble of molecules at a high temperature, which are moving around freely. Since physical systems tend towards lower energy states, the molecules are likely to move to the positions that lower the energy of the ensemble as a whole, as the system cools down. However molecules actually move to positions which increase the energy of the system with a probability 𝑒^{−Δ𝐸/𝑇′}, where Δ𝐸 is the increase in the energy of the system and 𝑇^′ is the current temperature. If the ensemble is allowed to cool down slowly, it will eventually promote a regular crystal, which is the optimal state rather than flawed solid, the poor local minima.

In function optimization, a similar process can be defined. This process can be formulated as the problem of finding a solution, among a potentially very large number of solutions, with minimum cost. By considering the cost function of the proposed system as the free energy and the possible solutions as the physical states, a solution method was introduced by Kirkpatrick in the field of optimization based on a simulation of the physical annealing process. This method is called Simulated Annealing. The Simulated Annealing algorithm to solve such problems is given below.(1)Start with some state, S.(2)𝑇^′=𝑇^′₀(3)Repeat {(4)While (not at equilibrium){(5) Perturb S to get a new state S𝑛(6) Δ𝐸=𝐸(𝑆_𝑛)-𝐸(𝑆)(7) If Δ𝐸<0(8)  Replace S with S𝑛(9) Else with probability 𝑒^{−Δ𝐸/𝑇′}(10)   Replace S with S𝑛(11) }(12) 𝑇^′=𝐶∗𝑇^′/∗ 0<C< 1 ∗/(13)} Until (frozen)

In this algorithm, the state, 𝑆, becomes the state (approximate solution) of the problem in question rather than the ensemble of molecules. Energy, 𝐸, corresponds to the quality of 𝑆 and is determined by a cost function used to assign a value to the state and temperature, 𝑇^′ is a control parameter used to guide the process of finding a low cost state where 𝑇^′₀ is the initial value of 𝑇^′ and 𝐶 (0<𝐶<1) is a constant used to decrease the value of 𝑇^′.

Procedures of SA Functions
(a)Representation. A “𝐾-dimensional real vector” S=(s1, s2,… s𝐾) is used to represent a solution, where 𝑠₁, 𝑠₂,…, 𝑠_𝐾 represent different decision variables of the problem under optimization.(b)Initialization. In this step an initial solution from the search space is generated. Different components 𝑠₁, 𝑠₂,…, 𝑠_𝐾 are randomly generated from their bounds such that constraints of the problem are satisfied. This solution is taken as initial state S in the above algorithm.(c)Perturbation Function. To find a neighbour solution of 𝑆=(𝑠₁,𝑠₂,…,𝑠_𝐾), a random integer 𝑟 is selected in the range [1…K]. Then 𝑠_𝑟 is replaced by randomly generated value within the boundary of the 𝑟th component of S or a random number 𝑟₁ is generated between −0.25 and +0.25 using random number generator, and 𝑠_𝑟 is replaced by 𝑠_𝑟+𝑟₁, so that S satisfies constraints of the problem. Final solution is taken as neighbour solution of S.(d)Energy Function. Value of the objective function (𝑓) due to solution 𝑆, 𝑓(𝑆), is taken as energy of 𝑆 if the problem is of minimization type otherwise −𝑓(𝑆) is taken as the energy function of the solution 𝑆.(e)Cooling Schedule. Initial temperature 𝑇^′₀ is taken according to different parameter values of the energy function, and reducing factor for 𝑇^′ (temperature), 𝐶 is taken as 0.999.