New Adaptive Surrogate-Based Approach Combined Swarm Optimizer Assisted Less Tuning Cost of Dynamic Production-Inventory Control System

Today, most manufacturing control systems are complex and expensive, so they are limited to employ a small number of function evaluations for optimal design. Yet, looking for optimization methods with the less-computational cost is an open issue in engineering control systems. This paper aims to propose an effective adaptive optimization approach by integrating Kriging surrogate and Particle Swarm Optimization (PSO). In this method, a novel iterative adaptive approach is utilized using two sets of training samples including initial training and adaptive sample points. The initial training points are designed by space-filling design, while the adaptive points are generated using a new jackknife resampling approach. The proposed approach can effectively convergence towards the global optimal point using a small number of function evaluations. The efficiency and applicability of the proposed algorithm are evaluated using the optimal design of the fractional-order PID (FOPID) controller for some benchmark transfer functions. Then, the introduced approach is applied for tuning the parameters and the sensitivity analysis of the FOPID controller for a dynamic production-inventory control system. The results are in good agreement with the results reported in the literature, while the proposed approach is executed with a lower computational burden.


I. INTRODUCTION A. STATEMENT OF PROBLEM AND MOTIVATIONS
Control engineering refers to the use of automatic control to make systems or processes reach the desired behavior while operating under certain constraints. This discipline has been being intensively enlarged over the past decades due to the advancement of modern technologies and the development of new systems, in particular intelligent systems [1]. Proportional-Integral-Derivative (PID) controllers with feedback control structures are widely used for industrial process The associate editor coordinating the review of this manuscript and approving it for publication was Sudipta Roy .
control. While more powerful control techniques are readily available, the PID controller is still popular due to its relative simplicity and applicability to a wide range of industrial control problems, see [2]- [5]. PID controller employs proportional, integral, and derivative gain parameters to reduce an error signal (e) to be closer to zero. The output of a PID controller equal to the control input of the plant in the time-domain is: where u (t) is the overall control function and K p , K i , and K d denote the proportional, integral, and derivative gain parameters, respectively. e (t) calculates an error value as the difference between the desired setpoint and a measured process output in the time moment t. The controller attempts to minimize the error over time by adjusting a control variable u (t). The appropriate function of PID strongly depends on the allocated value for each gain of K p , K i , and K d . Currently, fractional-order controllers are being extensively used by many scientists to achieve the most robust performance of the systems [6], [7]. The main reason for choosing Fractional-Order PID (FOPID) controllers is their additional degrees of freedom that result in a better control performance [2], [8]. A generalized FOPID controller was first introduced by [9] which proposed PI λ D µ controller-involving a λ order integer and a µ order differentiator. The differential equation of a fractional order PI λ D µ controller is defined by: The reliability of the FOPID controller depends on the optimal design of three gain parameters K i , K p , and K d and two order parameters (λ, µ). There are a few software tools available for designing and tuning fractional-order controllers. Tuning a fractional PID controller is a challenging task because there are five parameters to tune. This means that there are two more parameters as compared to a classical PID controller [2], [7].
The main motivation of the current study is to develop a new method for tuning FOPID controllers when the proposed method is specified in two main specifications including i) a black-box method that does need applying dynamic mathematical expressions of a control system, and ii) a less-expensive method with the lower number of function evaluations. Here, it should also be recalled the ''No Free Lunch'' (NFL) theorem presented by [10]. This theorem has logically proved that there is no best-suited optimization method for solving all optimization problems. In other words, a particular optimization method may show very promising results on a set of problems, but the same algorithm may show poor performance on a different set of problems [11]. This theorem also motivates the current study to develop a new optimization method for tuning the FOPID controller.

B. RELATED WORKS
Different tuning methods such as the Ziegler-Nichols method, auto-tuning, self-tuning method, and optimization-based tuning for the fractional-order controller have been proposed in the literature [1], [12], [13]. However, there is still a challenging research problem since most of the integer and fractional-order PID controllers are suffering from badly-tuned parameters, being the reason why they operate improperly in non-optimal regions [14], [15]. The optimization-based tuning methods for the complex control systems are classified into mathematical methods, simulation-based methods, and surrogates [16], [17]. Investigating less computationally expensive (i.e., a smaller number of iterations or function evaluations) methods for achieving optimal design controller has also become a main challenging topic in the control engineering practice. To overcome such computational difficulties, researchers have applied surrogate-based learning methods (e.g. polynomial regression, Kriging, and radial basis function) [18], [19]. The Kriging surrogate has been used as a widespread global approximation technique that is applied widely in control design problems [20], [21].
The best technique to adopt depends on the complexity level of the control system. Some other well-known metaheuristic methods applied in intelligent control are Grey Wolf Optimizer (GWO) [22], GA [23], [24], Particle Swarm Optimization (PSO) [24], [25], Ant-Colony Optimization (ACO) [26], [27], and Evolutionary Programming (EP) [24], Ant Lion Optimizer (ALO) [28]. The simplicity and global characteristics of such metaheuristics have been the main reasons for their extensive applications in off-line optimum control system design [29]. There are also some main modifications in the performance of mentioned optimizers applied in the engineering design, see [30]- [33].
Optimal control theory is a branch of mathematics developed to find optimal ways to control a dynamic system. One of the common control problems in engineering design applications is a production/inventory control system [34]. A block diagram representation of classical inventory and order-based production control system is schematically shown in Figure 1 [35], [36]. In the production-inventory control system, the controller aims to generate sophisticated decisions that control the order rate and the inventory level based on trade-offs between production-inventory cost and customer satisfaction [37], [38]. Order policy can be defined as PID controller systems. This kind of control for the production-inventory system is frequently used in a variety of applications [37], [39]. The application of the control theory in inventory-production systems has been illustrated in different studies, see [40]- [43].

C. MAIN CONTRIBUTIONS
The main contributions of our study can be summarized as follows: • In this paper, we introduce a new easy-toimplement algorithm for the optimal design of the VOLUME 9, 2021 FOPID controller. This optimization approach can adaptively and effectively converge to a global optimum point. In this approach, we employ the advantages of both the Kriging surrogate and PSO metaheuristic, combined with a new framework of the jackknife leaveone-out approach.
• The proposed method is inexpensive in terms of the required number of function evaluations for the optimization procedure. This method is proper for such expensive FOPID control systems when the model is limited to a small number of function evaluations (also called simulation experiments).
• A new dynamic model for a production-inventory control system is introduced. We also define a new objective function applied in FOPID optimization. The proposed optimization algorithm is used to tune the FOPID controller for such a production-inventory control system. Computational results indicate that the proposed algorithm outperforms alternative approaches from accuracy, robustness, and performance perspectives. Notably, as aforementioned, there are a variety of optimizers and modifications in the literature which have been developed to employ efficiently in different engineering design applications. However, the current paper limits its scope to compare the proposed algorithm with the five most common optimizers which frequently used in the tuning of integer and fractions PID controllers namely GA, ACO, PSO, GWO, and ALO. The rest of this paper is organized as follows. Section 2 provides the preliminaries required for our proposed algorithm. The algorithmic framework of the proposed approach is given in Section 3. The design procedure of the fractional PID controller for three typical benchmark transfer functions and a new framework for a dynamic production-inventory control system are presented in Section 4. Finally, this paper ends by presenting the limitations and advantages of the current study. Finally, some concluding remarks are mentioned in Section 5.

II. PRELIMINARIES
In what follows, first the basic definition of Kriging surrogate is given. Subsequently, the PSO method as a common optimizer in the control engineering practice is presented.

A. KRIGING
Since a Kriging was first introduced by Krige [44] in 1951, it has been used as a widespread global approximation technique that is applied widely in engineering problems [45]. Kriging is an interpolation method that can cover deterministic data and is highly flexible due to its ability to employ a various range of correlation functions. The higher accuracy of Kriging models comparing to the other alternatives such as response surface modeling is confirmed in the literature [46], [47]. In a Kriging model, a combination of a polynomial model and realization of a stationary point is assumed by the form of: The polynomial terms of f j (X ) are typically first or second order of response surface approach and coefficientsβ j are regression parameters (j = 0, 1, . . . , k). The term ε describes approximation error and the term Z (X ) represents the realization of a stochastic process which is the most time normally distributed Gaussian random process with zero mean, and variance σ 2 , and non-zero covariance. The correlation function of Z (X ) is defined by: where σ 2 is the process variance and R(x i , x p ) is the correlation function that can be chosen from different functions proposed in the literature. Due to the tuning of the correlation function with sample data, the Kriging is extremely flexible to capture the nonlinear treatment of the model. Among the studies in the literature, some studies have been found which sufficiently describe Kriging metamodel methodology, see [45], [48].

B. PARTICLE SWARM OPTIMIZER
The canonical PSO algorithm was proposed by [49] and is inspired by the social behavior of swarms such as bird flocking or fish schooling. The parameters of PSO consist of the number of particles, position of the agent in the solution space, velocity, and neighborhood of particles (communication of topology). The PSO algorithm begins with initializing the population. The second step is to calculate the fitness values of each particle, followed by updating individual and global bests as the third step. Then, velocity and the position of the particles become updated (step four). The second to fourth steps are repeated until the termination condition is satisfied [50], [51]. The velocity and position of the PSO algorithm are formulated as follows [49], [50]: where w is the inertia weight factor and v t id and x t id are particle velocity and particle position respectively. The expression d is the dimension in the search space, i is the particle index, and t is the iteration number. Expressions c 1 and c 2 represent the speeds of regulating the length when flying towards the most optimal particles of the whole swarm and the most optimal individual particle. The term p i is the best position achieved by particle i so far and p g is the best position found by the neighbors of particle i. The expression rand (0, 1) shows the random values between 0 and 1. The exploration happens if one or two of the differences are large, i) differences between the particle's best p t id and previous particle's position x t id , ii) differences between the population's all-time best p t gd and previous particle's position x t id . Besides, exploitation occurs when these two values are both small. PSO has attracted wide attention in control engineering design problems due to its algorithmic simplicity and powerful search performance [52]. However, a PSO algorithm that requires a large number of fitness evaluations before locating the global optimum is often prevented from being applied to computationally expensive real-world problems [53]. Therefore, surrogate-assisted PSO metaheuristic optimization algorithms have been focused on in the literature, see [21], [54].

C. CONVENTIONAL JACKKNIFE RESAMPLING
Jackknife was first introduced by Quenouille [55] and was named by Tukey [56]. The application of the jackknife method involves a leave-one-out strategy for the estimation of a parameter (e.g., the variance) in a dataset. In general, the jackknife resampling approach has been used to estimate surrogate prediction error and investigates new sample points to improve the accuracy of a surrogate [57]. This method is based on the leave-one-out approach, and it uses an existing set of data and is not required to re-run the expensive simulation model. Assume the initial set of training points is designed with n training points (i = 1, 2, . . . , n).
To select a new sample point between c candidate points (j = 1, 2, . . . , c), and add it into the set of training points, drop one point from the current set of training points, −i, and construct a surrogate based on n − 1 remaining points. Note that to avoid extrapolating by a surrogate (e.g., Kriging, radial basis function, or neural networks), the sample points on the vertices are not dropped. The jackknife's pseudo-value for candidate j is calculated as below: . . , n and j = 1, 2, . . . , c (7) whereŷ −0 j is the original prediction for candidate j with surrogate over the whole training sample points and y −i j is the prediction for candidate j with surrogate over n − i sample points (remove i th sample point from n set of points). The jackknife's variance s 2 j is computed for candidate j by employing relevant pseudo-values: So, the candidate with maximum s 2 j is selected as a winner and is entered in the set of current training sample points after computing its relevant true response with an original simulation. All the steps are repeated till stopping creation is satisfied. Among the previous studies, no specific stopping criterion was found to be appealed for in almost all cases [58]. However, it can be defined based on a limitation of computational time or cost.

III. PROPOSED ALGORITHM
In this section, the proposed promising optimization algorithm using hybrid Kriging surrogate and PSO metaheuristic based on adaptive improvement using the jackknife leaveone-out technique is explained. For this purpose, we first explain the objective function used in the proposed optimization model. Then, we sketch the algorithmic steps in the proposed approach. This approach is specified for the optimal design of the FOPID controller in complex and expensive systems.

A. OBJECTIVE FUNCTION DEFINITION
For an optimal design of a FOPID controller, a suitable objective function that represents system requirements should be defined based on some desired specifications. Some conventional output specifications in the time-domain are the overshoot (M p ), rise time (T r ), settling time (T ss ), and steady-state error (E ss ). It is worth noticing that using different performance indices makes different solutions for the optimal FOPID controllers [7], [52]. Time-domain performance criteria including IAE, ISE, and ITSE formulas are defined as follows: These performance criteria have advantages and disadvantages. For example, one disadvantage of IAE and ISE criteria is that their minimization can result in a response with a relatively small overshoot but a long settling time [52]. Although, the ITSE performance criterion can overcome the disadvantage of the ISE and IAE criteria [59]. Here, we propose and use a new cost function using all time-domain specifications including overshoot, rise time, settling time, and steady-state error besides the ITSE performance index as below: (12) where K = K p , K i , K d , λ, µ is the parameter vector and α∈ [0, 1] is the weighting factor. The optimum selection of α depends on the designer's requirements and the characteristics of the plant under control. To reduce the rise time and settling time, α can be set to be small (close to 0). Alternatively, if α is to be set larger (close to 1), then the overshoot and steady-state error are reduced more. In contrast, if α is set to 0.5, then all the performance criteria (i.e., the overshoot, rise time, settling time, and steady-state error) will have the same importance for optimization in the objective function of the FOPID model.

B. A NEW JACKKNIFE RESAMPLING CRITERION
The conventional jackknife approach in Eq.(7) considers the prediction error to select the best points among some candidate points and adding them adaptively into the set of training sample points to increase the accuracy of approximation by the surrogate. However, the conventional formulation in Eq. (7) investigates the candidate points with a higher jackknife variance (Eq.(8)), when it doesn't consider the possibility that the selected point between candidate points is located VOLUME 9, 2021 surrounding the global optimal point in a design space. In other words, the winner point (with higher jackknife's variance) between all the candidate points, may or may not have a lower fitness function in an optimization model under study. In this study, we aim to adaptively estimate the possible location and investigate the global optimal point in the design space inspired by the conventional formulation in Eq. (7). Accordingly, we propose a new formulation of the jackknife leave-one-out approach. Assume that the initial set of training points is designed with n training points (i = 1, 2, . . . , n); to select a new sample point between c candidate points (j = 1, 2, . . . , c), we apply a new jackknife leave-one-out formula as below: . . , n and j = 1, 2, . . . , c (13) where y −0 min is the smaller value of fitness function among all the points in the set of the current training points. This value is obtained using the output of the true simulation model for each sample point in the set of training points. The expression y −0 j denotes the prediction of the jth candidate point using the surrogate that is constructed overall training points. The expression y −i j shows the prediction of the jth candidate point using the surrogate that is constructed over the set of training points except for ith sample point. Between all j = 1, 2, . . . , c candidate points, the point with higher value L j is selected as a winner and is added to the current set of training points. Typically, Eq.(13) investigates the point in design space with a higher possibility of improvement compared to the other candidate points. Notably, since Eq.(13) employs a surrogate to predict output for each candidate point instead of running the true simulation model, it doesn't increase the computational cost and the required number of function evaluations for the optimization procedure. Accordingly, we can evaluate a larger set of candidate points in each iteration without concerning the computational cost. In this study, we employ the PSO optimizer to produce candidate points and investigate a point with maximum L j in Eq.(13).

C. ALGORITHMIC FRAMEWORK
In this study, we propose a new adaptive optimization algorithm by integrating the Kriging surrogate, PSO optimizer, and a new jackknife leave-one-out approach. The flowchart of the design procedure of the proposed algorithm is shown in Figure 2. Notably, this algorithm is low-cost to implement for such expensive systems when the system is limited to a small number of function evaluations required for the optimization model. Here, we specify the proposed algorithm to apply in a less expensive optimal design of the FOPID controller. In the FOPID optimization model, five decision variables including three gains K i , K p , K d and two order λ, µ parameters are considered. In this model, Eq.(12) is considered as the objective function. Algorithm 1 sketches the algorithmic framework and the main steps of the proposed approach.
In this algorithm, the set of training sample points for fitting Kriging and interpolating the whole design space includes i) initial training points and ii) update points. Here, we use a common space-filling design method namely Latin hypercube sampling (LHS) to design the initial training sample points. LHS was first introduced by McKay and colleagues [60]. It is a strategy to generate random sample points while guaranteeing that all the portions of the design space are depicted. In general, for n input variables, m sample points are produced randomly in m intervals or scenarios (with equal probability). We apply a PSO metaheuristic to investigate the whole design space to select the best sample point and update the training set which is applied to reconstruct the Kriging surrogate accordingly. In the proposed algorithm, to select the update point in each sequential run, the PSO optimizer looks for a point with a larger jackknife value in Eq. (13). These new update points can adaptively improve the accuracy of Kriging interpolation. This surrogate can interpolate the whole design space particularly locations with the higher possibility to surround the optimum points. For a typical transfer function as an instance, the initial train sample points, and adaptive sample points for the FOPID tuning problem are shown in Figure 3 using the proposed algorithm.
Notably, the optimization procedure for this typical transfer function including other instances will be discussed more in Section 4. As can be seen in this example, two sets of sample points are produced through the optimization procedure including initial sample points (e.g., here the initial sample size is designed with 15 points) and the adaptive sample points that sequentially are designed and added into the set of train points to update Kriging (here also 15 adaptive points are obtained through 15 sequential iterations). The optimum point obtained for this example is shown for each combination of decision variables K i , K p , K d , λ, µ . It should be noted that in this figure, the size of bubbles depicts the value of the objective function. It can be seen that the adaptive Algorithm 1 The Proposed Algorithm in Pseudocode Input: The upper and lower bound of decision variables (i.e. five decision variables for FOPID tuning is K = K p , K i , K d , λ, µ ), the parameter α ∈ [0, 1] in Eq. (12), and a total number of sequential iterations to update the model (h). Output: The optimal point found by the proposed approach. begin Step 1. Design initial training sample points using space-filling design methods (e.g., Latin hypercube sampling) with size n, i = 1, 2, . . . , n.
Step 2. Run the original model (i.e. here we use simulation model) for each sample point and obtain the relevant true output y i regarding the objective function (here C α (K ) as formulated in Eq. (12) is used as the objective function for the FOPID tuning problem).
Step 4. Fit a Kriging surrogate over sets of input/output data pairs. Step 5. Run the PSO optimizer and select the optimizer's result as a winner point: • Set the maximization of a jackknife leaveone-out criterion in Eq.(13) as the fitness function for PSO.
• Employ the constructed Kriging to obtain the predicted output (an approximation for the objective function of the optimization model) y j , (j = 1, 2, . . . , c) applied in Eq.(13). Step 6. Run the original simulation model and obtain the true simulation output for the selected winner sample point.
Step 7. Update the set of input/output data pairs (set of training sample points).
Step 8. Update Kriging surrogate over the new set of input/output data pairs. Step 9. Update h = h + 1.
Step 9. Repeat Step 5 till Step 9 until reach the total number of sequential iterations, h (stopping criterion).
Step 10 Select the point with minimum objective function (C α ) between all training sample points (including initial samples and h adaptive sample points) as the model's optimal point. end sample points with proper accuracy (lower objective function) can converge in the place that the final optimal result is located. In other words, the proposed algorithm can estimate the location of the optimum point with a small number of sample points (e.g., 30 sample points in this example) including initial sample points and adaptive (sequential) sample points. In the following section, we evaluate the applicability and effectiveness of the proposed algorithm for the optimal design of the FOPID controller for three typical transfer functions and one practical case in a dynamic production-inventory control system.

IV. NUMERICAL EXAMPLES A. OPTIMIZATION PROCEDURE SETUP
In this section, to verify the performance of the proposed algorithm, comparative experiments are carried out in the less-expensive (with the small number of function evaluations) optimal design of the FOPID controller for three typical transfer functions and a new dynamic structure of production-inventory control system. The proposed algorithm in this study is developed for such complex systems when the optimization model is restricted to employ a small number of function evaluations. To evaluate the performance and accuracy of the proposed algorithm in comparison with common intelligent control methods, we select five common optimizers that are widely used in the optimization of control systems including PSO [49], GWO [11], ALO [28], ACO [26], and GA [61]. For designing an optimal FOPID controller, determining the decision variables K = K i , K p , K d , λ, µ associated to the minimization of the objective function C α (K ) (see Eq. (12)) is the main issue. Here, the proposed optimization algorithm is performed to minimize the objective function C α (K ). For this purpose, step response of the plant is used for computing all performance components in Eq.(12) including overshoot (M p ), rise time (T r ), settling time (T ss ) and steady state error (E ss ), and ITSE. Here, we assume that the optimization model is limited to only 30 function evaluations (simulation experiments). Moreover, we only let the proposed algorithm and other five optimizers PSO, GWO, ALO, ACO, and GA to employ only 30 function evaluations to obtain the optimal result for optimal design of FOPID controller. We simulate all control plants in Matlab R /Simulink environment. The size of the time-step is fixed at 0.01 and the time-domain (simulation time) is considered T = 20. The setpoint for all the problems is fixed to 1. Simulink does not have a library for the FOPID. Therefore, the controller from the library of FOMCON: a Matlab R toolbox for fractional-order system identification and control [62] which allows for the computation of the fractional-order derivative and integration is used. It is assumed that the suitable ranges of the control parameters for all control problems are as follows: and λ ∈ [0, 1] .
It should be noted that the parameter α in Eq.(12) is set to α = 0.5 for all control plants (in this case, all the performance criteria have the same merit in the objective function).
To study and compare the effect of randomness in the proposed algorithm with five common optimizers including PSO, GWO, ALO, ACO, and GA, we repeat each stochastic optimization method 10 times for each problem and compute the statistical measures including mean, Standard Deviation (Std), and the Signal to Noise Ratio (SNR) as the robustness criterion as below: where µ and σ show the mean and Std of optimizer obtained by repeating the optimization model (here 10 separate repetitions). Since we performed a minimization of the model's output (see the objective function in Eq. (12)) to find the optimal parameters of the FOPID controller, the formulation of SNR in Eq. (14) has the opposite sign by Taguchi formulation. Additionally, a weighting parameter ω is introduced to allow for individual emphasis on the minimization of variations. The smallest value of SNR in Eq.(14) depicts the better stochastic optimization method with higher accuracy (smaller mean of simulation output) and higher robustness (smaller standard deviation) against randomness in the model.

B. THREE TYPICAL TRANSFER FUNCTIONS
Three different control plants are extracted from [52] with transfer functions as below: The following process is performed to determine the optimal design of the FOPID parameters (K i , K p , K d , λ, and µ) using the proposed algorithm for the mentioned control plants. First, we design a small set of sample points (here 15 sample points) using the Latin Hypercube Sampling (LHS) method (here we use the Maltlab R function ''lhsdesign'' to construct LHS design). Then, we run the original simulation model to obtain relevant outputs (see Eq. (12)). Next, we fit a Kriging surrogate over a set of input/output data pairs (train sample points). After that, we run the PSO optimizer to select another sample point to add to the set of train sample points. The maximization of L j criterion in Eq.(13) is set as the fitness function for the PSO optimizer. The constructed Kriging is used to estimate the model's output used in Eq. (13). The result gained by PSO is selected as a winner point. We run the original simulation model again to obtain the true output in the winner point. The input-output data pair for the winner point is added to the existing set of training sample points. We also update the Kriging surrogate accordingly. This adaptive iteration is repeated until the stopping criterion is satisfied. In this study, we stop adaptive sampling when the model reaches 15 iterations (and accordingly 15 adaptive sample points are designed). Therefore, a total of 30 function evaluations regarding the designed sample points (initial and adaptive points) are performed through the implementation of the proposed algorithm for the optimal design of the FOPID controller. To make a fair comparison between the proposed algorithm and the other common FOPID solvers, we only let each optimizer e.g., PSO, GWO, ALO, ACO, and GA evaluate the main objective function 30 times. For this purpose, we set 3 initial populations and 10 maximum iterations for PSO, GWO, and ALO optimizers, 3 ants with 10 iterations for ACO, and a population size of 5 with 6 generations for GA. Other parameters are adjusted as follows. Min and max inertia weights, velocity clamping factor, cognitive constant (c1), and social constant (c2) for PSO is fixed to 0.7, 0.9, 1.5, 2, and 2 respectively. For the ACO optimizer, the evaporation rate and the number of nodes for each parameter are adjusted to 0.7 and 10,000 respectively. The obtained results for the plants G 1 (s), G 2 (s) and G 3 (s) are summarized in Tables 1, 2 and 3, respectively. The pairwise comparison with statistical measures including mean, Std, and SNR, and the unit step responses of the control plants for the control plants G 1 (s), G 2 (s) and G 3 (s) are shown in Figures 4-9, respectively. It is seen that the proposed algorithm to tune the FOPID controllers outperform those of the PSO, GWO, ALO, ACO, and GA obtained controllers. In all the obtained results in figures, it can be observed that the proposed algorithm gives higher accuracy (lower cost function C α (K ), see Eq.(12)) comparing to the other five methods. Besides, the proposed algorithm shows more robustness comparing to the other stochastic optimizers with a smaller standard deviation for results obtained over 10 repetitions of the method. From the simulation results, it can be found that the proposed algorithm applied in the tuning of FOPID controllers produces the smooth curve for the output in conjunction with little fluctuation and small overshoot.

C. PRODUCTION-INVENTORY CONTROL SYSTEM 1) DYNAMIC OF PRODUCTION-INVENTORY CONTROL SYSTEM
A good implementation of the control theory applied to production-inventory control has been presented by [36], [63]- [65]. Here, we consider a factory producing single homogeneous goods and having a finished goods warehouse. TABLE 1. FOPID optimal results in G 1 (s) using different methods for 10 repetitions (α = 0.5).
In this study inspired by [35], we design the dynamic of the production-inventory control model as below: (18) where I (t) indicates the inventory level at time t, and δ(t), (0 ≤ δ ≤ 1) shows a portion of the inventory that cannot be used in the next time step due to some reasons like damage or frustration. The optimal order policy O (t) is designed by the FOPID controller, see Figure 10. The expression D(t) shows the demand on time step t. In this model, we consider the failure rate of the product that is increased by the passing of time. So, we assume that the production rate follows a nonlinear exponential distribution, β = exp(−θ.t), (0 ≤ θ ≤ 1). It means that in each time step, the production process only produces the portion of β, (0 ≤ β ≤ 1) from the order policy.   Thus, the production rate P (t) at time step t is computed by P (t) = β.O(t). The block diagram representation of the proposed dynamic production-inventory control system is shown in Figure 11.

2) RESULTS AND DISCUSSION
To design an optimal FOPID controller based on the proposed algorithm, we implement the simulation model using Eq. (18) in the Matlab R , Simulink environment. In the current instance, we assume that the demand rate follows a nonlinear exponential distribution as D (t) = 2exp(0.05t). Besides, the failure rate is assumed θ = 0.02, so production rate equals to β = exp(−0.02.t). We also assume that in the whole simulation period, the constant rate of 10 percent (δ (t) = 0.1) from remaining inventory in time step t cannot  be used in the next time step t + 1 due to damage or frustration. The optimization process for the production-inventory control plant is done the same as what was done for three typical plants in the previous section. To determine the optimal values of the FOPID gain and order parameters (K = K i , K p , K d , λ, µ ) by the proposed algorithm, we first design 15 sample points using the LHS method. Each sample including the vector of K = K i , K p , K d , λ, µ is sent to the Simulink block and the values of five performance criteria in the time-domain (M p , T r , T ss , E ss , and ITSE) are calculated. Accordingly, the objective function C 0.5 (K ) (see Eq.(12) with α = 0.5) is computed for each sample regarding the obtained performance criteria.
Afterward, a Kriging surrogate is fitted over the set of input/out data pairs (15 samples). We employ the PSO optimizer to design a new sample point and update the set of input/output data pairs accordingly. For this purpose, we set the maximization of the proposed leave-one-out jackknife criteria regarding Eq.(13) as the fitness function in the PSO optimizer. The constructed Kriging interpolation is applied to estimate the value of the objective function required to FIGURE 11. Block diagram representation of the FOPID controller in the dynamic production-inventory control system. compute L j (Eq.(13)) as the fitness function in PSO. As we employ Kriging instead of the original simulation model to compute the objective function, we are not worried about the computational cost with PSO optimizer in our algorithm (thus a large number of the initial population and the maximum number of iterations can be considered for the swarm optimizer). The PSO used in the proposed algorithm has the same parameter adjustment as the PSO that individually is used for direct FOPID tuning. Finally, the point with maximum L j is selected as a winner and this new adaptive sample point is added to the set of input/output data pairs which is employed to update the Kriging accordingly. This adaptive sampling process is continued until the algorithm reaches the stopping criterion. Regarding the proposed algorithm (please see the pseudocode in ''Algorithm 1''), we stop the algorithm (terminate repeating Step 5 till Step 9) when the total number of sequential iterations (h) is reached the predefined number (stopping criterion). Here for the current test functions, we stop the algorithm when 15 adaptive points are designed (thus 15 sequential iterations are derived). Finally, among all the 30 sample points (15 initial samples and 15 adaptive sample points), the best point with a lower objective function C 0.5 (K ) is returned as the global best solution for K . For comparison, we also solve the FOPID tuning problem in production-inventory control system with five common optimizers of PSO, GWO, ALO, ACO, and GA. Note that to provide a fair comparison, we also let each optimizer employ 30 function evaluations to search for the optimal result. The parameter adjustment for each optimizer is the same as what was mentioned in the previous subsection for tuning three typical transfer functions. The comparative results for 10 individual repetitions of each method are summarized in Table 4.
Besides, Figure 12 shows statistical results including the boxplot, mean, standard deviation, and SNR values over the obtained results using each method.
It is seen that the FOPID tuning based on the proposed algorithm outperforms those of the FOPID tuning based on other optimizers from the accuracy (lower objective function C 0.5 (K )) point of view. Besides, the proposed algorithm outperforms other optimizers in terms of robustness   (smaller standard deviation) against randomness regarding different repetitions of stochastic optimization methods. On the other hand, as can be seen from Figure 12, the proposed algorithm has a smaller SNR measure for different values of parameter ω in Eq. (14) comparing to the other optimizers. Figure 13 reveals the step response of the VOLUME 9, 2021 TABLE 5. The statistical results of sensitivity analysis for the parameter α in Eq. (12) in production-inventory control system. The results obtained over 10 repetitions using the proposed algorithm. motor rotation for different FOPID controllers. Obviously, the FOPID controller tuned by the proposed algorithm produces a smooth and fast output comparing to the controllers tuned by the PSO, GWO, ALO, ACO, and GA methods.

3) SENSITIVITY ANALYSIS
Finally, the sensitivity analysis is performed with different values for the parameter α in Eq. (12) to investigate the effects of the parameter α on the optimal results obtained by the proposed algorithm in the dynamic production-inventory system. To do this, the parameter α is changed from 0 to 1 with a step size of 0.25. The obtained results are given in Table 5. The statistical results for 10 repetitions of the proposed algorithm regarding each value of α are revealed in Figure 14.
As can be seen from the results, the cost function C α (K ) is reduced by increasing α from zero to one. It also can be seen that by increasing α from zero to one, the overshoot (M p ) and steady-state error (E ss ) are reduced and rise time (T r ) and settling time (T ss ) are increased. In general, the controller tuned by α = 0.75 has a smaller overshoot and steady-state error, in contrast with the case α = 0 which has a lower rise time and settling time. In addition, the best and worst ITSEs are obtained in the cases of α = 0.75 and α = 0.25 respectively.

V. CONCLUSION
In this paper, a new simple modified hybrid approach is proposed that integrates Kriging surrogate and particle swarm optimization to overcome the computational cost (high number of function evaluations) required in the common optimization methods in optimal design of FOPID controller. We propose a new criterion for adaptive sampling strategy inspired by a leave-one-out jackknife point of view. In the proposed algorithm, we consider increasing the accuracy of the surrogate in regions of design space that the global optimal point is located. The competitive results between the proposed algorithm and with the five most common optimizers in the tuning of control systems namely PSO, GWO, ALO, ACO, and GA show that the proposed method can provide a good balance between exploration and exploitation to the avoidance of local optima. The advantages of this methodology are its less computation burden as well as high-quality solution and speedy convergence. The results also revealed that the proposed algorithm is more robust in the stochastic environment of the optimization model. This paper has some limitations including: • The proposed method employs a Kriging surrogate to train the model using space-filling sampling strategies. Therefore, the approximate errors cannot be ignored when solving simulation-based optimization problems particularly with complex functions and nonlinear structures. A challenge for optimization under restricted budgets will be to find • the right degree of approximation (smoothing factor) from a relatively limited number of samples [66].
• In this study, the application of the proposed algorithm is evaluated and compared with other optimizers for the optimal design of the FOPID controller in three benchmarks and the dynamic production-inventory control plants. However, the performance of the proposed algorithm can be examined for more challenging cases with a higher level of complexity and nonlinearity.
• Other new optimizers in the literature with modified formulation and structure compared to the five classic optimizers applied in this study for comparison of results can be employed to evaluate more the efficiency of the proposed hybrid surrogate-metaheuristic method in the tuning of the controller with a small number of function evaluations. Therefore, future research will be devoted to overcoming the aforementioned limitations of the current study.