A New Simulation-Based Robust Optimal Design of FOPID Controller for Five Bar Linkage Robot in a Cyber-Physical System

: This paper aims to further increase the reliability of optimal results by setting the simulation conditions to be as close as possible to the real or actual operation to create a Cyber-Physical System (CPS) view for the installation of the Fractional-Order PID (FOPID) controller. For this purpose, we consider two different sources of variability in such a CPS control model. The first source refers to changeability of a target of control model (multiple setpoints) because of environmental noise factors and the second source refers to an anomaly in sensors that is raised in a feedback loop. We develop a new approach to optimize two objective functions under uncertainty including signal energy control and response error control while obtaining the robustness among the source of variability with the lowest computational cost. A new hybrid surrogate-metaheuristic approach is developed using Particle Swarm Optimization (PSO) to update the Gaussian Process (GP) surrogate for a sequential improvement of the robust optimal result. The application of efficient global optimization is extended to estimate surrogate prediction error with less computational cost using a jackknife leave-one-out estimator . This paper examines the challenges of such a robust multi-objective optimization for FOPID control of a five-bar linkage robot manipulator. The results show the applicability and effectiveness of our proposed method in obtaining robustness and reliability in a CPS control system by tackling required computational efforts. Multiple Setpoint.


Introduction
Nowadays, developing processes in the engineering world are strongly associated with computer simulations. These computer codes can collect appropriate information about the characteristics of engineering problems before actually running the process.
Computer simulations can provide a rapid investigation of various alternative designs to decrease the required time to improve the system. In addition, most numerical analyses for engineering problems make a well-suited use of mathematical programming. The main goals of simulation include what-if study of a model or sensitivity analysis and optimization and validation of the model [1]. The essential benefit of simulation is its ability to cover complex processes, either deterministic or random, while eliminating mathematical sophistication [2]. Clearly, because of the complexity of mathematical formulation analyzing in many real-world optimization problems, simulation-optimization methods become necessary to find more interest and popularity than other optimization methods [3]- [5]. CPS combine physical objects or systems with integrated computational facilities and data storage [6]. CPS is a key enabling technology in systems intelligence. In CPS embedded computers and networks, the physical processes are controlled usually with feedback loops where physical processes affect computations and vice versa [7]- [10].
CPS are multidimensional and complex systems that integrate the cyber world with the dynamic physical world. Integrating physical processes with computer systems is the main challenge presented in CPS as the computational cyber part continuously senses the state of the physical system and applies decisions and actions for its control [11].
In industrial practice, many CPS systems have been designed by decoupling the control system design. In this way, CPS and real-time interaction are achieved in order to monitor and control physical entities in a reliable, safe, collaborative, robust, and efficient way [12], [13]. Using precise calculations to control a seemingly unpredictable physical environment is a great challenge [31]. After the CPS control system is designed and modeled by extensive simulation, tuning methods need to be expanded to address uncertainty and random disturbances in the system. In addition, ignoring the impact of uncertainty on the optimization model, the obtained optimal results may be far from the true optimum settings [32]. One of the main features in a reliable CPS design is stability feature (robustness), which means no matter how the environment generates noise and uncertain factors, the control system should always reach a stable decision

Control
Cyber Space result eventually [33]. Robustness in CPS control system seeks to achieve a certain level of performance with possible modeling errors in the forms of parametric or nonparametric uncertainties [34]. However, considering uncertainty and random disturbances, while keeping the function and operation of the system, has been computationally time-consuming and costly.
Because of uncertainty, more complexity in real-time control implementation of CPS is unavoidable. So, looking for less expensive computational methods of optimization considering uncertainties has become interesting among most engineering applications.
Surrogate-based methods can 'learn' the problem behaviors and approximate the function value. These approximation models can accelerate the function evaluation as well as the estimation of the function value with an acceptable accuracy. In addition, they can improve the optimization performance and provide a better final solution.
In this paper, a new outline of robust real-time optimization in the CPS control model under the effect of environmental factors (also known as noise factors or uncertainty, see [4], [5]), and variability in feedback loop due to sensor's anomaly is studied. The main contributions of this study are as following: i. We propose a new hybrid surrogate/metaheuristic method combining the GP surrogate and PSO algorithm. We apply the PSO metaheuristic to update the GP surrogate for sequential investigation of a robust optimal result. The proposed hybrid GP/PSO algorithm has the advantages of both GP surrogates in learning the behavior of the model in an efficient global optimization with PSO metaheuristic in convergence searching for optimum results.
ii. We apply the straightforward jackknife leave-one-out technique to estimate surrogate prediction error applied in efficient global optimization. This method can estimate surrogate prediction error apart from the type of surrogate using a training set of sample points.
iii. An augmented bootstrapping technique is used to analyze the sensitivity of robust optimization results. This technique can apply the initial set of Input/Output (I/O) data instead of resampling the model.
iv. An exhaustive search method for optimization under uncertainty with a straightforward procedure is developed. It can be easily applied in real operation of the CPS framework. In this study, the proposed algorithm is applied to provide robust tuning of the FOPID controller over two sources of variability.
The first source is related to real-time setpoint that is predicted by learning from collected data (e.g. surrogate) over CPS environmental factors and the second variability is found in output's feedback due to anomaly in sensors.

v.
Energy consumption and response error are optimized as a robust multiobjective optimization model by Pareto frontier estimation in the real-time computational part of the CPS model (see Figure 1).
The rest of this paper is organized as follows. Section 2 provides more details about realtime FOPID control when two types of uncertainties (noises) including environmental factors and sensor anomaly are considered in a CPS framework. Materials and methods of the proposed algorithm to handle robust multi-objective optimization of a CPS control system are elaborated in Section 3. In Section 4, the applicability and effectiveness of the proposed approach are examined to provide robustness and reliability in the robust optimal design of the FOPID controller in the CPS framework of a five-bar linkage robot manipulator. Finally, this paper is concluded in Section 5.

Problem Statement
The existing uncertainties and anomalies in the cyber environment have resulted in emerging concerns about the traditional control system [34]. In real-time control of CPS, physical process variables are monitored and processed by intelligent controllers for keeping the values of safety parameters between the given thresholds.
Environmental conditions can affect the system dynamics and also the controller function [9]. The precision of computing must interface with the uncertainty and the noise in the physical environment [44]. The physical world, however, is not entirely predictable. It is normal that the CPS do not operate in a controlled environment. So, it must be robust to uncertainty (unexpected conditions) and adaptable to subsystem failures [8].

FOPID controller
In this paper, for better control, fractional-order controller is used. Currently, fractional-order controllers are being extensively used by many scientists in order to achieve the most robust performance of the systems [45]. The main reason for choosing FOPID controllers is their additional degrees of freedom that result in a better control performance [46], [47]. A generalized FOPID controller was first introduced by [48] which proposed controller involving a order integer and a order differentiator. The differential equation of a fractional-order controller is defined by: The reliability of FOPID controller depends on the optimal design of three gain parameters ( , , ) and two order parameters ( , ). However, we try to further increase the reliability of the tuning result by setting the traits of the simulation model to be as close as possible to the practical condition to make a CPS outline for the FOPID controller. The FOPID control system with a single setpoint does not express the aspects of the behavior that are essential to the system in the context of CPS. Moreover, we challenge the robust control to achieve CPS stability when the uncertainty in environmental conditions is the source of variability of the setpoints in the control system. In addition, an uncertain anomaly in sensors causes the noise (variability) in the control feedback loop. Moreover, we aim to tune the FOPID controller robustly in such a CPS control system with real-time setpoints and noise in the model's feedback. Figure 2 shows the control outline of CPS with real-time setpoints and noise in the model's feedback. The application of the integer-order and fractional-order of the PID controller in CPSs has been studied in [47], [49]- [51].   . Anomaly in sensor's feedback is function of uncertain variable ̃. FOPID gain parameters and order parameters are tuned robustly in such a way to make CPS insensitive against sources of variability in system.

Objective functions
This study aims at optimizing a robust multi-objective model of the FOPID tuning in the CPS framework by considering two different objective functions (e.g. performance criteria). The first objective function is targeted to manage the Signal Energy Control (SEC) that is consumed in the time domain 0 ≤ ≤ as follows: and where 1 is a user-defined big value that is used for normalizing the first objective function in [0,1], so that 1 > log ( max 0≤ ≤ ).
We define the second objective function based on the Response Error Control (REC) (i.e. inspired integrated absolute error) as below: ). Notably, we use a logarithmic scale for both objective functions to smooth the large differences between the values (i.e. cases in which one or a few points are much larger than the bulk of the data). As mentioned earlier, the real-time setpoint ( ) in Eq.(5) can be predicted on-time by an easy-to-apply supervised learning like polynomial regression as a function of environmental uncertain factor(s).

Overall objective function
In order to combine both objective functions including the signal energy control (see Eq. (2)) with response error control (see Eq.(4)) to be used as a single objective model, we apply -mertic approach by = 2 (i.e. for more information about -mertic approach in multi-objective optimization, refer to [52]). Assume = (1,2, … , ) is the vector of input combination, then we define the Overall Function (OF) as below: where is a user-defined weighting factor (0 ≤ ≤ 1) that indicates the tendency of the model toward optimization based on each objective function 1 and 2 , see Eq. (2) and Eq.(4). Fluctuating this magnitude ( ) provides the capture of Pareto frontier (also called Pareto optimal efficiency) to make a trade-off between each objective function.
This approach is a classical method to solve optimization problems when the model is faced with multiple criteria [53]. In fact, the set of optimal solutions obtained from fluctuating in [0,1] provides an estimate of the Pareto frontier.

Proposed Algorithm
In this section, we propose a promising technique for optimization under uncertainty using augmented efficient global optimization using the jackknife leave-one-out technique to estimate GP prediction error hybrid GP/PSO method. For this purpose, we first explain the main materials and methods used in the proposed algorithm briefly and then sketch the algorithmic steps in the proposed approach.

Gaussian process (GP) surrogate
GP, that is also known as kriging, is a non-parametric Bayesian approach to supervised learning [54]. GP is an interpolation method that can cover deterministic data and is highly flexible due to its ability to employ various ranges of correlation functions [55].
In a GP model, a combination of a polynomial model and the realization of a stationary point are assumed by the form of: where the polynomial terms of ( ) are typically the first or the second-order response surface approach and coefficients ̂ are regression parameters ( = 0,1, … , ). This type of GP approximation is called the universal GP, while in the ordinary GP, instead of ( ), the constant mean = ( ( )) is used. The term describes the approximation error and the term ( ) represents the realization of a stochastic process which in general is a normally distributed Gaussian random process with zero mean, variance 2 , and non-zero covariance. The correlation function of ( ) is defined by: where 2 is the process variance and ( , ) is the correlation function that can be chosen from different correlation functions which were proposed in the literature (e.g. exponential, Gaussian, linear, spherical, cubic, and spline), see [56], [57]. Today, GP surrogate has been used as a widespread global approximation technique that is applied widely in control engineering design problems [38], [58].

Particle swarm optimizer (PSO)
The canonical PSO algorithm was proposed by [59] and is inspired by the social 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 [60], [61]. The PSO algorithm is formulated as follows [59]- [61]: where and are particle velocity and particle position respectively. is the dimension in the search space, is the particle index, and is the iteration number.
Expressions 1 and 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 is the best position achieved by particle so far and is the best position found by the neighbors of particle . The expression (0,1) shows the random values between 0 and 1. The exploration happens if either or both of the differences between the particle's best ( ) and previous particle's position ( ), and between population's all-time best ( ) and previous particle's position ( ) are large. In addition, the 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 [62], [63]. However, 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 [64]. Therefore, surrogate-assisted PSO metaheuristic optimization algorithms have been focused in the literature, see [64]- [66].

Uncertainty management
Here, we follow [37], [67], [68] and inspire Taguchi's overview of robust design [69] for dealing with uncertainty as a source of variability in the model. However, we expand Taguchi robust design terminology and apply its definition for environmental noise factors in such a CPS control system. But in this study, we replace statistical approach of Taguchi viewpoint with augmented efficient global optimization using jackknife leave-one-out technique and hybrid GP/POS approach. Furthermore, we first intersect two sampling design sets. One sampling design is for decision variables (inner array) and another is for uncertain variables (outer array). Given that = (1,2, … , ) is the vector of sample points over decision variables, and = (1,2, … , ) is the vector of uncertainty scenarios, so × input combinations are designed, and the real model (or true simulation model) are evaluated × times to collect relevant simulation outputs, see Figure 4. Assume is the × matrix of simulation outputs (i.e. in this study the simulation outputs include OF values that can be obtained regarding Eq. (6)), thus mean and standard deviation (Std) for each arrow in can be computed by the following equations: Signal-to-Noise Ratio (SNR) as introduced by Taguchi [69], [70] is a robustness criterion based on the mean and the Std of a system response . Given that, is the smaller the better type, Taguchi assumed zero as the minimal possible response value. Accordingly, he formulated the following SNR as the robustness criterion: Since we performed a minimization of the model's output (here is overall function, see Eq. (6)) to find the optimal parameters of FOPID controller, the formulation of the SNR in Eq.

Efficient global optimization using a jackknife leave-one-out strategy
A common formulation of efficient global optimization has been developed in the outline of the expected improvement criterion, see [71], [72]. The expected improvement (EI) method has been developed in engineering design problems to adaptively improve the local and global search of optimal points (i.e. control a trade-off between exploration and exploitation properties) [73]. This method has been combined with two main parts. The first statistical part consists of the design of experiments and surrogate techniques and the second part involves evolutionary algorithms. If is considered for the arbitrary point , an improvement function over the best point that is so far computed with is defined as formulation of efficient global optimization in term of expected improvement creation is constructed as below:  shows the value of the best signal to noise ratio that is obtained from true data of the original simulation model, and ̂ is GP surrogate prediction on candidate point . The terms Φ and ∅ depict the cumulative distribution function (CDF) and probability density function (PDF) of a standard normal distribution respectively.
The first phrase (Φ) in Eq. (14) is related to the local search and the second phrase (∅) is related to a global search. In the search for the next best point among all the candidate points, the point with maximum EI term in Eq. (14) is selected and replaced with the best point so far obtained. This procedure is continued until − 0 ≤ , where is a user-defined threshold, or an allocated computational cost (e.g. fixed number of repetitions) is reached. However, to find the neighbor points (candidate points) around the current best point, different strategies of sampling design methods can be used such as factorial designs [74] and space-filling design [75]. Here, we apply PSO global optimizer to investigate the maximum EI among the whole design space.
In order to estimate the surrogate prediction error for th candidate point (̂) in Eq. (14), simulation experiments can be resampled [71], [72]. The authors in [75] [77]. 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 [78]. In this study, we are motivated to use the jackknife leave-one-out approach to estimate surrogate prediction error (̂) required in Eq. (14) formulation, because this method uses an existing data and does not require to re-run the expensive simulation model. Here, this method is used to predict GP prediction error while it can be used for other surrogates as well. Let ̂ denotes the prediction of GP surrogate that fitted over all samples (input combinations), therefore the GP perdition error in th candidate point (̂) can be estimated through the jackknife leave-one-out approach as the steps in Algorithm 1.
Input: Set of input combinations and relevant output (SNR). Output: Estimation of surrogate prediction error for th candidate point. begin Step 1. Select samples from the complete set of combinations ( = 1,2, … , ) when = − and is a set of samples located in vertices (i.e. we aim to avoid extrapolating of GP surrogate).
Step 4: Predict output for th candidate point (̂− ) using the GP surrogate constructed from the previous step.
Step 5: Implement three previous steps for all samples computing relevant predictions.
Step 6: Apply the jackknife estimator to obtain the estimation of surrogate prediction error for th candidate point as below: In this study, we develop a new hybrid surrogate/metaheuristic method applied in robust efficient global optimization and optimization under uncertainty. We apply a PSO metaheuristic to update a GP surrogate for sequential investigation of a robust optimal point. The proposed algorithm can handle robust efficient global optimization by the exhaustive search method that can be applied in real operation of CPS control frameworks. The algorithmic representation of the proposed approach is presented in  Cross two DOE sets over decision and uncertain variables using space filling design, ( × ) combinations.
Fit a GP surrogates, over SNR values of overall function.
Run PSO for maximization of EI criterion (or minimization of -EI) and obtain winner point.
Run CPS control model in the winner point for uncertainty scenarios and collect relevant outputs for each objective function.
Compute overall functions in the winner point.
Compute mean, SD and SNR of overall function for the winner point.
Update the GP surrogates of SNR by adding the input/output data for the winner point on the set of training sample points.
Obtain robust optimal point.
By varying in [0,1], the estimation of Pareto frontier is obtained.

Finish
Algorithm 2 Proposed robust simulation-optimization approach.
Output: The estimation of the Pareto frontier by a set of robust optimal points found by the proposed approach. begin Step 1. Design crossed array (using the space-filling design) by crossing two sets of experiments with dimensions × as below: • An inner array matrix with dimension (l × n x ) where l is the number of sample points for decision variables and n x is number of decision variables (e.g. in FOPID tuning n x =5 decision variables including three gain K i , K p , K d and two order λ, μ parameters).
• An outer array matrix with dimension (m × n z ) where m is the number of sample points (uncertainty scenarios) for n z uncertain variables (e.g. here in represented CPS control system n z = 2 including ŝ(t) and α).
Step 5. Fit a GP surrogate over sets of I/O data (with input combinations and relevant values).
Step 6. Define an initial best point among the set of I/O data obtained from Step 4 (the point with the smallest SNR regarding Eq. (13)).
Step 7. Set expected improvement criterion (see Eq. (14)) as an objective function in PSO optimizer algorithm (i.e. with minimizing of − ( )) and obtaining a winner point.
Step 8. Run the real CPS model (e.g. original simulation model) in the winner point for combinations of uncertainty (scenarios) designed in Step 1 and obtain relevant outputs for each objective function.
Step 9. Obtain OF values for the winner point regarding uncertainty scenarios.
Step 10. Compute mean and Std of the winner point using Eq.(11) and Eq.(12).
Step 11. Update the set of I/O data = (1,2, … , + ), when is the number of the sequential runs.
Step 12. Fit a new GP surrogate over an updated set of I/O data (with + training points and SNR as outputs).
Step 13. Update (if needed) the best point obtained so far to a point with smallest SNR ratio among all the sample points (including initial training points and points which are added so far for updating of surrogate, see Step 11) and repeat Step 7 till Step 12 until stopping rules are satisfied (e.g. stop sequential updating if − 0 ≤ , or ≤ , where and are user-defined thresholds).
Step 14. If stopping rule(s) is satisfied, then set the best point obtained so far as a robust optimal point of the model. The best point so far has the smallest SNR value among all sample points including initial samples and updating sample points.
Step randomly in intervals or scenarios (with equal probability). Inspired by [80] in the case of non-independent multivariate input variables, the desired correlation matrix can be used to produce distribution-free sample points in LHS. For more information, refer to [37], [81]. In the represented CPS control system in this study, outputs for two

Sensitivity analysis of optimal results using augmented bootstrapping approach
In this study, the main idea behind the proposed algorithm is to perform sensitivity analysis to expand the information obtained from robust efficient global optimization.
Estimating a single optimal point using a particular response may be inaccurate because of variability in the surrogate. Thus, we derive a series of possible responses that take into account a degree of uncertainty by providing confidence regions or prediction intervals. The author in [82] has mentioned two alternative strategies for bootstrapped resampling as follows: • In each set of bootstrapping, both sets of input (design) combination ( ) and noise (uncertain) combination ( ) are resampled randomly.
• The resampling is adapted to noise or uncertain component ( ) only while keeping the deterministic input combination ( ) fixed.
Here, to find the bootstrapped set of data, a model is resampled times ( = 1,2, … , ) The superscript '*' is a common symbol for bootstrapped values [56]. The expression gives two-sided CIs. Bonferroni's inequality suggests that Type I error rate for each interval per output is divided by the number of outputs (here is SNR). If the values of bootstrap estimate ( + ) * are sorted from low to high, then ⌊. ⌋ and ⌈. ⌉ respectively denote floor and ceiling function to achieve the integer part and round upwards.
In this study, as inspired by [68], [84], the particular augmented bootstrapping approach is used for costly simulation running. In such a case, assume the set of sample points is fixed and only old data to fit surrogate with enough replication is available and new simulation replicating is very expensive. This augmented bootstrapping approach does not imply extra computational cost because of resampling and required simulation running to find a bootstrapped set of data. ( = 1,2, … , ) denotes the set of sample points and each is repeated times ( = 1,2, … , ). We assume that the original set of data obtained from the original simulation model is available (size × ) when is the number of scenarios for uncertainty and is the number of input combinations.
Moreover, the augmented bootstrapping procedure is sketched in Algorithm 3.
Algorithm 3 The augmented bootstrapping procedure.
Input: Set of I/O data, and robust optimal point. Output: Estimation of CIs. begin Step 1. Set = 1 and = 1.
Step 3. Replace the th original output , (selected from the old data) with the bootstrap output , * = , * .
Step 8: Compute bootstrapped CIs using Eq.(15) for the robust optimal point obtained by the proposed algorithm as elucidated in Section 3.2). end Note that, regarding the Step 1 till Step 5, it can be seen that a random number with replacement in [1, ] is selected and regarding the selected number, we choose the relevant response in an outer array (see the structure of the crossed array design explained in Section 3.1.3) that was previously collected from original simulation model and has the same column number. For the same input combination, we repeat this procedure times and collect different responses or may have the same responses (i.e. because the random selection is done with replacement). This procedure is also repeated for other input combinations. Therefore, the data matrix with row and column is constructed.

Numerical Example
Here, the proposed algorithm is specified for the robust optimal design of FOPID controller in CPS control of five-bar linkage robot manipulators. In the continue, we first explain the dynamics of the five-bar linkage robot manipulators. Next, the robust optimal design of FOPID controller in the CPS framework of a five-bar linkage robot manipulator is obtained using the proposed algorithm in this paper.

Dynamics of five-bar linkage robot manipulator
It should be noted that 1 depends only on 1 but not on 2 . Similarly, 2 depends only on 2 but not on 1 . This discussion helps to explain the popularity of the parallelogram configuration in industrial robots. If 3 2 3 = 4 1 4 , then two angles 1 and 2 can be adjusted independently without worrying about interactions between the two angles.

Robust FOPID tuning in CPS control of five-bar linkage robot manipulator
Here, the main goal is to obtain a robust optimal design of FOPID controller in such a CPS control model as elucidated in Section 2. We simulate the five-bar linkage robot manipulator using Eq. (16) in Matlab®/Simulink environment. Simulink does not have a library for the FOPID. Therefore, the controller from the library of FOMCON: a Matlab® toolbox for fractional-order system identification and control [89], which allows for the computation of the fractional-order derivative and integration is used.
Numeric values of the parameters of the five-bar manipulator dynamics are taken from [62], [90] as shown in Table 1. From the data driven in Table 1, it is revealed that 3 2 3 = 4 1 4 , thus we can perform robust optimal design of controller for one motor, where the same results are also valid for the second motor. In FOPID tuning, five decision variables including , , , , and are considered as decision variables.
The search procedure for the robust optimal result is done in the ranges as [62]:  Two performance criteria are considered as outputs of the model including Eq.(2) and Eq. (4)  Moreover, according to an updated set of training samples, a new GP surrogate is constructed after each sequential EI iteration. It is important to note that we avoid extrapolation of GP surrogate in each sequential iteration by setting two different rules, i) we consider a death penalty for any point that is investigated by PSO and is located out of bounds of training points, ii) in order to estimate GP prediction error using jackknife leave-one-out approach, we only remove input combination's rows that don't locate on the margin of design space (see Section 3.1.5). The obtained results on performing the proposed algorithm are discussed in the following section.

Results and discussion
We perform the proposed algorithm for three different values of = 0.25, = 0.5, and = 0.75 in computing OF (see Eq. (6)). In addition, to evaluate the effect of randomness in sampling design methods, each optimization set was repeated 5 times. Table 2, Table   3, and Table 4  update combinations that are crossed by 9 uncertainty scenarios). We consider all the three best points over all five repetitions as robust optimal points for the FOPID controller using the proposed algorithm for each = 0.25, = 0.5, and = 0.75 separately. Figure 7 shows the magnitudes of the EI criterion and the best SNR obtained by sequential expected improvement over five different repetitions of the proposed algorithm for = 0.25, = 0.5, and = 0.75. In addition, mean and Std of OF related to the best point so far (smaller SNR) in each sequential EI iteration are shown in Figure   8. It should be considered that two stopping rules are adjusted. The rules mention that EI value becomes smaller than 0.01, or reaches 15 sequential iterations. Figure 9 includes   Table 3 Robust FOPID optimal results using proposed algorithm over five different repetitions for = 0.5.  Table 4 Robust FOPID optimal results using proposed algorithm over five different repetitions for = 0.75.  Then, we compare the results obtained by a proposed algorithm with three common FOPID optimization methods in the literature ( [47], [62], [92], [93]) including the PSO metaheuristic [59] (i.e. when directly used in optimization procedure), Grey Wolf

Repeat
Optimizer (GWO) [94] and Ant Lion Optimizer (ALO) [95]. To make a fair comparison, we run true simulation model regarding each uncertainty scenario (total 225 simulation runs for each set of FOPID optimal point). Table 5, Table 6, and Table 7 include the results of comparison using proposed algorithm and three global optimizers. As can be seen, the proposed algorithm can obtain smaller SNR (more robustness) comparing to others. In addition, the proposed algorithm obtains smaller mean and Std of overall functions for = 0.25, = 0.5, and = 0.75 when it shows the effectiveness of the proposed algorithm in making a tradeoff between both REC and SEC objectives. Figure   9 indicates the variability of the model in each optimal point among simulation outputs including OF, REC, and SEC. Variability in the CPS control model is observed because of uncertainty (source of variability). In Figure 10, the variability of the model is evaluated over 225 different combinations (scenarios) of both uncertain factors including ( ) and ̃.   In order to analyze the sensitivity of robust optimal results obtained by the proposed algorithm and estimate the variability which occurred due to randomness in designing sample points, we used the augmented bootstrapping method explained in Section 3.3.
Here, based on the obtained results from original GP surrogate, the FOPID parameters    bootstrapped GP surrogate and 95% confidence intervals (CIs) over robust optimal point obtained by original GP surrogate for = 0.5. Augmented parametric bootstrapping is performed using on hand set of input/output data provided among original optimization program. bootstrapped GP surrogate and 95% confidence intervals (CIs) over robust optimal point obtained by original GP surrogate for = 0.75. Augmented parametric bootstrapping is performed using on hand set of input/output data provided among original optimization program.

Conclusion
In this paper, a new hybrid surrogate/metaheuristic-based robust simulationoptimization algorithm is proposed that possesses the advantages of both GP surrogate in learning the behavior of the model in efficient global optimization and PSO metaheuristic in convergence searching of optimum results. We smooth the application of PSO using GP surrogate instead of the original simulation model to diminish computational cost due to a large number of fitness evaluations required for the global optimizer when used individually. Also, this simple modified algorithm is developed in such a way to handle computational complexity to obtain optimal and robust FOPID design in the CPS control system. In such a system, we also consider conflict between multiple objective functions and uncertainty in the model. Here, we apply this approach to the robust optimal control design of a five-bar linkage robot manipulator to depict the applicability and effectiveness of the proposed algorithm. Comparative simulation results reveal that the proposed hybrid GP/PSO-based robust efficient global optimization algorithm can effectively and robustly tune the parameters of the FOPID controllers. From an application point of view, the introduced technique is simple and fast and has a suitable control on error and energy of system and it can be easily implemented in real-world applications of CPS control systems.