Robust Optimization with Interval Uncertainties Using Hybrid State Transition Algorithm

: Robust optimization is concerned with ﬁnding an optimal solution that is insensitive to uncertainties and has been widely used in solving real-world optimization problems. However, most robust optimization methods suffer from high computational costs and poor convergence. To alleviate the above problems, an improved robust optimization algorithm is proposed. First, to reduce the computational cost, the second-order Taylor series surrogate model is used to approximate the robustness indices. Second, to strengthen the convergence, the state transition algorithm is studied to explore the whole search space for candidate solutions, while sequential quadratic programming is adopted to exploit the local area. Third, to balance the robustness and optimality of candidate solutions, a preference-based selection mechanism is investigated which effectively determines the promising solution. The proposed robust optimization method is applied to obtain the optimal solutions of seven examples that are subject to decision variables and parameter uncertainties. Comparative studies with other robust optimization algorithms (robust genetic algorithm, Kriging metamodel-assisted robust optimization method, etc.) show that the proposed method can obtain accurate and robust solutions with less computational cost.


Introduction
Real-world optimization problems are subject to uncertainties due to, for example, the presence of uncontrolled changes in environmental conditions [1,2], the lack of complete knowledge of models [3,4], and the manufacturing tolerances on actual processes [5]. According to the classification of optimization problems under uncertainties [6], deriving an optimal solution that is insensitive to uncertainties is defined as robust optimization (RO).
An optimal solution is robust if, when uncertainties exist, the values of the corresponding objective functions and constraint functions fluctuate within acceptable ranges [7]. By adjusting the decision variables to counteract the effects of uncertainties, the performance of optimal solution can be improved and the degradations of objective or constraint function values can also be avoided. In general, there are two methods to describe the uncertain parameters in the optimization problem: probabilistic models and nonprobabilistic models. Probabilistic uncertainties are customarily based on the statistic information, such as the mean and variance, and they are usually handled by optimizing the expected value of the solution [8]. Moreover, probabilistic uncertainties can also be modeled by other methods, such as fitness approximation [9], multiobjective approach [10,11], and sampling-based methods [12]. Since it is difficult to obtain the accurate probability distribution information of uncertain parameters and the optimal solution can hardly guarantee complete robustness, the applications of probabilistic robust optimization methods have been limited. Nonprobabilistic uncertainties methods are usually based on interval uncertainty modeling [13,14], evidence theory [15], and possibility theory [16]. In this paper, the uncertainties are modeled as a certain interval that can be obtained without a presumed probability analyzes the robust performance of the proposed method based on a comparative study of eight number of algorithms on seven examples. Section 5 concludes this paper and discusses the directions of future research.

Robust Optimization Problem
In general, the deterministic optimization problems can be defined as follows: s.t. g i (x, p) ≤ 0, i = 1, · · · , n, where f (·) is the objective function and g(·) is the constraints function, and n is the number of constraints. The vector x is the decision variable whose lower and upper bounds are x l and x u , respectively, and p is the parameter of the problem. In optimization problems, uncertainties can be involved in both decision variables and parameters. Thus, the formulation of the optimization problem under interval uncertainty is given as where [x] and [p] are interval numbers corresponding to the uncertain decision variables and uncertain parameters. They can be expressed as where x c and p c are the nominal value of x and p, respectively, with ∆x and ∆x being the lower and upper bounds of decision variable (x) variation, and ∆p and ∆p being the lower and upper bounds of parameter (p) variation. For simplicity, it is usually assumed that the nominal value is the central value of the variation range, which implies that ∆x = −∆x and ∆p = −∆p.
For evaluating the solutions of the robust optimization problem in (2), three indexes are introduced:

1.
Objective robustness: This index, denoted as η f , is a measure of the sensitivity for the objective function to uncertainties. When decision variables and/or parameters fluctuate in their uncertain intervals, the objective function variations should still within an acceptable range. In engineering problems, the acceptable range of objective function is usually defined by decision makers according to design requirements.

2.
Feasibility robustness: This index, denoted as η g , is a measure of the sensitivity for the constraints to uncertainties. When decision variables and/or parameters fluctuate in their uncertain intervals, the constraints still should be satisfied.

3.
Optimality: This index, represented as f , is the objective function value. For a deterministic optimization problem, the optimum should be the solution with best objective value.
Based on the above indexes, the optimization problem in (2) can be reformulated as follows: g i (x c , p c ) ≤ 0, i = 1, · · · , n, where ∆ f 0 means the acceptable variation range of the objective function.
The above formulation about the robust optimization problem contains a nested double-loop optimization structure. The outer optimization functions to find promising nominal values of decision variables and the inner optimization is used to verify the robustness of candidate solutions. The nested double-loop optimization structure incurs high computational costs. Thus, this paper proposes a robust optimization method that solves the computationally costly optimization problem using (i) the state transition algorithm with sequential quadratic programming, called the hybrid state transition algorithm (H-STA), to solve the outer optimization problem, and (ii) the second-order Taylor series expansion to estimate the robustness indexes for inner optimization. In the sections that follow, the two techniques of the proposed RO method are further discussed.

Taylor Series Surrogate Model
In general, the inner optimization is performed iteratively for every candidate solution in the outer optimization. To reduce the computational cost, the procedure for inner optimization should be as simple as possible. In the proposed RO method, a second-order Taylor series surrogate model is used to approximate the objective function and calculate the extreme points for inner optimization.
Based on the Taylor's theorem, a multivariable function f (x, p) can be expanded by the Taylor series around (x, p) = (x c , p c ) Let then (5) can be written as The remainder term of the second-order Taylor expansion of f (x, p) around point (x c , p c ) can be written as follows: where θ ∈ (0, 1). It is worth noting that R(∆x, ∆p) is a function of a cubic polynomial with respect to ∆x and ∆p, and the values of ∆x and ∆p are usually small. Meanwhile, the higher-order derivatives have relatively small values compared to the lower-order derivatives [38]. The above two points make R(∆x, ∆p) a tiny value and guarantee the accuracy of the second-order Taylor expansion alternative model.
For solving the maximization problem in (4), the second-order Taylor series is adopted and the inner optimization problem can be transformed to where Similarly, the feasibility robustness index can be transformed to where In (9) and (11), the extreme points can be computed by the quadratic formula, and the maximum is calculated by backsubstituting into the original function in (4). With the calculation of the quadratic function, the computational cost of the inner optimization problem can be reduced.

State Transition Algorithm
For outer optimization problems, most metaheuristic methods have competitive performance. The state transition algorithm (STA) [34,39] is an intelligent optimization method based on the control theory of state space representation. The unified form of the generation of solutions in the STA method can be described as follows: where x k represents a state, corresponding to a candidate solution of the problem; u k is a function of historical states; A k and B k stand for state transition matrices; and y k means the fitness value of the objective function f . In STA method, there are four state transformation operators that generate candidate solutions: where (14)- (17) are the rotation transformation, translation transformation, expansion transformation, and axesion transformation, respectively. The parameters α, β, γ, and δ represent the transformation factors. The parameters R r , R t , R e , and R a are the random matrix with specific elements. The rotation transformation is a local search operator and the translation transformation has a function of line search. The expansion transformation is used for global search and the accession transformation is designed to strengthen singledimensional search ability. For a given solution, the aforementioned state transformation operators are performed alternately to generate candidate solutions. In general, these four operators in the STA method can find promising candidate solutions and converge to the global optimum point. However, in a robust optimization problem, the solution requires not only to be a superior candidate in the deterministic condition but also to satisfy the robustness requirements. Thus, it is important to strengthen the search ability and design appropriate selection mechanism when solving a robust optimization problem. Since the STA method is a stochastic algorithm which does not use the gradient information, the local search ability is restricted and the precision of the solution still needs some improvements. In this paper, a hybrid state transition algorithm that combines the improved STA operator with a traditional local search procedure is proposed to address the robust optimization problem.

Hybrid State Transition Algorithm for Robust Optimization Problem
In the STA method, the expansion transformation, as the main global search operator, still requires further improvement to enlarge the range of global search. In addition, the local search direction of the STA method is stochastic, which may have a slow convergence rate; thus, sequential quadratic programming (SQP) is used to exploit the local area and improve the precision of solutions. In this paper, the hybrid state transition algorithm (H-STA) for robust optimization is proposed and it combines the improved STA and SQP to maximize their advantages of global optimization and minimize their disadvantages of premature convergence. Moreover, in order to balance the feasibility, robustness, and optimality, an efficient selection mechanism is proposed to evaluate candidate solutions and select the best one as the final result.

Exploration Stage-Improved STA
The first stage of the H-STA method for robust optimization problems is to explore all search areas and find some promising candidates. In the basic state transition algorithm, there are four different transformation operators and the expansion transformation operator is the main global search operator. As shown in (16), the expansion transformation operator includes a Gaussian distribution matrix R e , which means that it can generate elements between [−∞, +∞] with probability. However, (16) also shows that the search range of expansion transformation not only depends on the expansion factor (γ) and the mean and standard deviation of R e , but also relates to x k . Thus, if the value of x k is small, the search range will be small. For example, if we set the initial point as [10,10] and [1,1] separately, and the lower bound and upper bound of x are set to −10 and 10, respectively, then the parameter setting of the expansion transformation operator is the same as in previous papers [39], which are γ = 1, and the mean and the standard deviation of R e equal 0 and 1, respectively. In this study, we use the expansion transformation operator to generate 500 candidates, and if the candidate value is out of the range, a random value within the range will be selected as a substitute.
The performance of the expansion transformation is shown in Figure 1, with Figure 1a showing the expansion transformation operator that can generate a candidate in the search space with the initial point [10,10]. However, when the initial point is set to [1,1], the search range of the expansion transformation operator becomes narrow (see Figure 1b). Thus, the global search ability of expansion transformation still requires further improvements.
One solution is to take into account the ranges of the decision variables in the expansion transformation: where R x = (x u − x l )/2 is the search radius of the decision variables, including the upper bound and lower bound of the decision variable. The pseudocode of the new expansion transformation operator is shown in Algorithm 1. We use the same parameter settings and initial points to perform the new expansion transformation operator, and the results are shown in Figure 2. It is shown that no matter how the initial point changes ( [10,10] in Figure 2a and [1,1] in Figure 2b), the candidates generated by the new expansion transformation can always distribute throughout the entire search space.
The structure of the improved STA method is shown in Algorithm 2. Firstly, the parameters are predefined and the initial solution is generated randomly. The parameter SE is the search enforcement, which means that every transformation operator will generate SE candidate solutions. During the optimization process, the rotation factor α decreases exponentially from a maximum value (α max ) to a minimum value (α min ) with the base fc. If the rotation factor α is less than α min , it will return to the maximum value α max .

Require:
oldBest: the best solution in the last transformation Ensure: Best * : optimal solution Using the improved STA (I-STA) method, the four transformation operators can explore the whole space for neighborhoods that may contain the global optimal solution. Once the STA method converges to a promising candidate and stops at this point after several iterations, then the SQP method is used to improve the precision of the solutions.

Exploitation Stage-SQP
The second stage of the H-STA method for robust optimization problems needs to exploit the local area based on a promising initial point. Sequential quadratic programming is effective for solving a nonlinearly constrained optimization problem since the goal of SQP is to find a locally optimal solution to the problem, which means that the algorithm will continuously optimize the local area near the current solution in order to reduce the value of the objective function as much as possible. That ensures the accuracy of SQP. As an iterative procedure, the SQP method transforms a nonlinear optimization problem into a quadratic programming subproblem, and by solving that subproblem, the solution will converge to a local minimum. Thus, with the promising candidate found by the STA method and taking it as the initial point, the SQP method will exploit the neighborhood around that point and can improve the accuracy of the solution. The SQP structure [40] is shown in Algorithm 3. Best k+1 = Best k + d k 8: k = k + 1 9: end if 10: Go to Step 1

Algorithm 3 Pseudocode of SQP method
We now explain the searching process of the SQP method; based on the Lagrangian function and quasi-Newton method, at each iteration, the outer optimization problem in (4) can be approximated by the following quadratic programming subproblem: where H k is the Hessian matrix of f (x) and d is the step length. For a given x k c , functions (19) can be calculated. In order to find a desired step length, the subproblem can be rewritten as By solving (20), the optimal step length is obtained. After the iterative calculation, the SQP method will stop when the optimal step length approaches 0.
Since the SQP method has a good exploitation ability and a low computation cost, it can effectively improve the accuracy of the solutions obtained by the I-STA method. Moreover, the I-STA method can generate a good initial solution for the SQP method, avoiding the possibility of falling into a local optimum due to a poor initial solution.

Selection Mechanism
Since the robust optimization problem requires to meet not only the constraints but also the robustness requirement, a selection mechanism is needed to balance the feasibility, robustness, and optimality of the solutions. To quantify the performance of the candidate solutions, two definitions are given as follows.

Definition 1 (Constraint violation). For the constraints g
For a solution x, if G(x, p) = 0, then x is a feasible solution; otherwise, x is an infeasible solution, and the larger the value, the stronger the constraint violation. The selection steps for two solutions of the robust optimization problem are as follows:
The feasible solution is always preferred to the infeasible solution; 2.
If two solutions are both infeasible, then the one having a smaller constraint violation value is preferred; 3.
If two solutions are both feasible, then the robustness indexes will be calculated by inner optimization, and the following criteria will be considered: The robust solution is always preferred to the nonrobust solution; (b) If two solutions are both nonrobust, then the one having a smaller robustness violation value is preferred; (c) If two solutions are both robust, then the one having better objective function value is preferred.
The code structure of the selection mechanism is shown in Algorithm 4. After generating candidate solutions in the outer optimization process, the inner optimization is performed only when the feasible solution exists. Note that the infeasible solution can never satisfy the requirement of feasible robustness. Best * : optimal solution

Framework of the H-STA for Solving Robust Optimization Problem
The framework of the H-STA method for robust optimization problem with interval uncertainty is shown in Figure 3.
Step 1 (Initialization): The first step includes the generation of the initial solution and initialization of the parameters.
Step 2 (Outer optimization -I-STA): Based on the initial solution, the I-STA method generates candidates using four transformation operators.
Step 3 (Selection mechanism): The proposed selection mechanism is used to select the best solution among many candidates. If there is a feasible solution in the candidate solutions, the inner optimization is conducted based on Taylor series surrogate model. Step 4 (Switching criterion): The switching criterion is used to determine whether the I-STA method should be replaced by the SQP method. The I-STA method is a stochastic optimization algorithm, and in the later iterations, the rate of convergence is slow. Nevertheless, the SQP method offers fast convergence rate leading to the optimal solution based on a good initial solution. Thus, if the objective value in the improved STA method changes very slowly, the SQP method will be performed to improve the rate of convergence. In this paper, the switching index is defined as follows: where f k and f k−1 are the objective function values of x k and x k−1 . If λ k is less than a threshold value (λ), the rate of convergence for the I-STA method is considered slow and the SQP method is carried out.
Step 5 (Outer optimization-SQP): Considering the final solution of the I-STA method as an initial point, the SQP method is performed to improve the precision of the solution.
Step 6 (Selection mechanism): The proposed selection mechanism is used to compare the final solution and the initial solution of the SQP method.
Step 7 (Stopping criterion): We use the maximum number of iteration as the stopping criterion. It is worth noting that the number of iterations in the SQP method is also taken into account when calculating the total number of iterations.

Verification Examples
In this section, the proposed method is applied to seven optimization examples with interval uncertainty. Table 1 gives the uncertainty occurrences in each example. To demonstrate the effectiveness of the H-STA method for robust optimization (H-STA-RO), the following methods are used for comparison: the I-STA method for robust optimization problems (I-STA-RO), the basic STA method for robust optimization problems (STA-RO), the SQP method for robust optimization problems (SQP-RO), and five well-known robust optimizers (the robust optimization method with Chebyshev surrogate models (I-RO) [22], the robust genetic algorithm (GA-RO), the Kriging metamodel-assisted robust optimization method (IK-GA-RO) [41], the robust optimization method based on Benders decomposition (BD-RO) [27], and the robust optimization method using differential evolution and sequential quadratic programming (DE-SQP-RO) [42]).

Uncertainty Examples 1, 2, and 3 Example 4 Examples 5 and 6 Examples 7 Occurrences
The complexity of H-STA-RO can be represented by the value of function evaluations (FE), which is calculated as where SE is the search enforcement; iter STA and N STA are the number of generations and transformation operators used in the I-STA method, respectively; the parameter FE SQP is the total function evaluation value in the SQP method. Since the inner optimization is only performed when the feasible solution is found, N R means the number of times that inner optimization performed, and R f and R g present the average value of function evaluations for candidates to obtained their objective robustness index and feasibility robustness index, respectively. In order to evaluate the performance of different methods, all results are obtained after 20 runs under MATLAB R2016a, Windows 10 machine with 2.40 GHz Intel core i5 and 16.0 GB RAM. The SQP method and the genetic algorithm method are performed by using "fmincon" and "ga" function, respectively, and all the parameters are set by the default values. The parameters included in the H-STA-RO method are selected empirically based on numerous experiments and application cases. In the standard continuous STA, the parameter settings are given as follows: α min = 10 −4 , β = 1, γ = 1, δ = 1, SE = 30, and fc = 2. Many numerical experiments and engineering applications have shown the effectiveness of the above parameter settings [34,35,39]. In this paper, to obtain better results for different problems, the parameters of α max , λ, and iter max are fine-tuned based on the following guidelines: • The rotation factor α, which controls the search range of the rotation transformation, is bounded as α min ≤ α ≤ α max . A larger value of α allows more explorations of the local search space, and a smaller value of α can refine the quality of solutions. The value of α max is typically set as 1 based on the previous study. However, for the problem in which the ranges of decision variables are less than 1, it is useless to search in a hypersphere with a radius equal to 1. Thus, the parameter α max in Example 4 (0 ≤ x 1 , x 2 , x 3 , x 4 ≤ 1) is adjusted according to the statistical analysis. As shown in Figure 4, by performing 20 trails in each test, we compare the average iterative results with different α max value. Given the same initial solution, the iterative curve with α max = 0.1 has fastest rate of convergence to search better solutions. Therefore, α max in Example 4 is set to 0.1.

•
The threshold value of switching index λ controls the frequency of the switching between two search operators. If one operator is trapped in a slow convergence, another search operator is taken into consideration. As shown in Figure 5, a larger value of λ can increase the switching frequency but it may give a low-quality solution under the SQP method (e.g., λ = 10 −2 ). A smaller value of λ may cause slow convergence (e.g., λ = 10 −5 ). In this paper, λ, as the threshold value of the relative difference between two objective values, is adjusted between [10 −4 , 10 −3 ].

•
The maximum number of iterations iter max depends on the complexity of the problem. For the two engineering optimization problems considered in this section, a choice of 80 to 100 iterations may be sufficient. For the four numerical problems, iter max is set to 40 to 60. Take Example 4 for example ( Figure 5), the red dotted line nearly has no update from 45 generation, thus iter max is set to 60 so that the convergence performance of the H-STA method is guaranteed.

Example 1
This unconstrained optimization problem is used to verify the accuracy of the inner optimization method. The uncertain problem is given by as follows: With uncertainties existing in the decision variables, the optimal solution should be the minimum of its objective upper bound, which can be represented as The upper bound of the objective function f u is calculated by inner optimization. Table 3 shows the best results obtained by three different robust optimization methods from 20 trails. The interval robust optimization (I-RO) method [22] uses the multi-island genetic algorithm (MIGA) as the outer optimization method, and its inner optimization is based on the second-order Chebyshev surrogate model. The linearization robust optimization (L-RO) process [22] replaces the inner optimization method of I-RO with the first-order Taylor series surrogate model. In the H-STA-RO method, the inner optimization method is based on the second-order Taylor series surrogate model. Table 3. Performance comparison of Example 1. To obtain the accurate value of the objective upper bound, we used the Monte Carlo method in the uncertain range around the design point (using 2 × 10 8 samples), which can be used as the reference value (R) for assessing the accuracy of inner optimization method.

I-RO (Second-Order Chebyshev) [22] L-RO (First-Order Taylor Series) [22] H-STA-RO (Second-Order Taylor Series)
The results show that the values of f u based on the second-order Taylor series model and the second-order Chebyshev model are closer to the reference value (R). Figure 6 illustrates the optimization results of these three methods. Figure 6a shows the contour lines of the objective function values, and Figure 6b is the results of the Monte Carlo test with 200 samples. We observe that the H-STA-RO method can find the global optimum and its inner optimization method provides sufficient accuracy to evaluate the robustness index.

Example 2
This example is a nonlinear numerical problem with uncertainty in decision variables [42]. The deterministic version of this example is a classical multimodal function called "peaks function". In the deterministic problem, the optimal solution is x = [0.2283, −1.6255] and the objective function value is f = −6.5511. When the decision variable x 1 is subject to interval uncertainty, the robust optimization problem takes the form: Based on the proposed selection mechanism, shown in Table 4 are the optimal results obtained by the DE-SQP-RO method, the GA-RO method, the STA-RO method, the I-STA-RO method, and the H-STA-RO method. Since this example is a classical multimodal function, it has a known local optimum and global optimum. To verify the global search ability of the I-STA method, the success rate p s of removing the local optimum (the percentage of successful runs in total runs) is analyzed.
In Table 4, the decision variable, the constraint function value, and the robustness violation value all correspond to the results with the best objective function value in 20 runs. The values for the success rate p s and the robust rate p r (the percentage of robust runs in total runs) are obtained by statistical analysis. The value for FE is the average value of the function evaluations and the standard deviation is also presented. The value of T represents the average runtime (in seconds). From Table 4, we observe the following: (1) The success rate of removing the local optimum in the STA-RO method is 60%, and there are five results that fall into the local optimum [−0.2606, 0.4667] with the objective function value 0.7881. In the I-STA-RO method, however, all the tests can find the results that are close to the global optimum. Thus, the new expansion transformation operator offers a better global search ability.
(2) The robust rates of the STA-RO method, the I-STA-RO method, and the H-STA-RO method are all 100% (with the proposed selection mechanism), whereas only 95% and 75% of the results in the GA-RO method and the SQP-RO method (without the proposed selection mechanism) are robust, which demonstrates that the proposed selection mechanism can obtain a robust solution with higher probability.
(3) In the DE-SQP-RO method [42], although its objective function value is better than that of the H-STA-RO method, the result cannot satisfy the requirements of the robustness according to the value of R.
(4) The average function evaluations and runtime of the H-STA-RO method are smaller than that of others (except the SQP-RO method); this is because on the one hand, the proposed selection mechanism can avoid useless calculations in inner optimization, and on the other hand, the SQP method can improve the efficiency of the search process. The proposed selection mechanism may also cause the variation of inner optimization computational cost, leading to a high standard deviation of FE.
To verify the robustness of the obtained solution in the H-STA-RO method, the Monte Carlo simulation is conducted. By using 200 samples around the nominal value within the uncertainty interval, the objective robustness index and the feasibility robustness index are calculated. In Figure 7, when decision variable x 1 is perturbed, the variations of the objective function (shown in Figure 7a) and the constraints (shown in Figure 7b) for the solution obtained by the H-STA-RO method are always within the acceptance range, whereas the deterministic solution cannot satisfy the robustness requirements.

Example 3
This constrained nonlinear problem originated from [41]. The problem formulation is:  Table 5. From Table 5, it is observed that the H-STA-RO method can find the optimal solution with 100% robustness, and its function evaluations and runtime are less than that of the GA-RO method, the STA-RO method, and the I-STA-RO method. Although the function calls of the IK-GA-RO method is smaller, the optimality and robustness of the optimum obtained by the IK-GA-RO method are inferior to that of the H-STA-RO method. Figure 8 shows the robustness indexes of the deterministic optimum and robust optimum obtained by the H-STA-RO method. The deterministic optimum violates the robust requirement at some points, but the robust optimum of the H-STA-RO method can remain stable within the uncertain range.

Example 4
This example [27,43] illustrates the solving of the robust optimization problem with uncertainty in both decision variables and parameters. The formulation of this problem is min Without incurring the uncertainties, the deterministic optimal solution is x = [0.5, 0.5, 0.5, 0.5] with f = 9.770. Table 6 shows the results for this example using robust optimization methods. These methods are the BD-RO method, the GA-RO method, the SQP-RO method, the STA-RO method, the I-STA-RO method, and the H-STA-RO method. Figure 9 shows the results of the Monte Carlo tests. From Table 6, it is observed that the objective function values of the I-STA-RO method and the STA-RO method are close to the optimal value, but it is hard to improve the precision. The SQP-RO method has a good rate of convergence but it cannot ensure the robustness of the solutions. Thus, the H-STA-RO method takes advantage of I-STA and SQP to obtain the robust solution and search the global optimum with less computational cost. For the GA-RO method, it has better objective function value but the computational efficiency and the robustness of the solution need to be improved. The BD-RO method can find the optimum with fewer function calls but its outer optimization approach is based on the gradient information; thus, its results are highly influenced by the initial point. Figure 9 also demonstrates that when there are uncertainties in both decision variable and parameters, the solution obtained by the H-STA-RO method can still satisfy the constraints.

Example 5: Welded Beam Design
The welded beam design ( Figure 10) is a classical constrained optimization problem in engineering applications [36,44]. The design objective is to minimize the total cost of the welded beam f . The decision variables are the thickness of the weld x 1 , the length of the welded joint x 2 , the width of the beam x 3 , and the thickness of the beam x 4 . The decision variables must satisfy the constraints about the shear stress (τ), the bending stress (σ), the buckling load on the bar (P c ), the end deflection of the beam (δ), and the side constraints. The deterministic optimal result is x = [0.2053, 3.2604, 9.0366, 0.2057] and f = 1.6956. When the decision variables are subject to uncertainties, the problem is modified as follows: 6 psi, E = 30 × 10 6 psi, P = 6000 lb, L = 14 in, τ max = 13600 psi Table 7 shows the results obtained by the DE-SQP-RO method [42], the GA-RO method, the SQP-RO method, the STA-RO method, the I-STA-RO method, and the H-STA-RO method. From Table 7, the H-STA-RO method can generate competitive solutions when compared with other solutions. This is because (1) based on the I-STA method, the H-STA method can search the whole space and choose promising candidates, and (2) the SQP method can search the local area and improve the precision of the solution. The selection strategy and the inner optimization method proposed in this paper can guarantee the robustness of the best solutions. The findings presented in Figure 11 verify the robustness of the final solutions with respect to their objective function ( Figure 11a) and constraints (Figure 11b).

Example 6: Pressure Vessel Design
The optimal design of pressure vessel (see Figure 12) is a typical application example to illustrate the min-max problem. To minimize the total cost of the vessel f , there are four decision variables to optimize: x 1 (the thickness of the shell), x 2 (the thickness of the head), x 3 (the inner radius), and x 4 (the length of the vessel without the head). The best deterministic result is f = 5886.4544, corresponding to x = [0.7785, 0.3848, 40.3389, 199.7753]. With the uncertainties in the decision variables taking into account, the pressure vessel optimization problem is formulated as follows: where x 1c and x 4c are the nominal value of x 1 and x 4 , respectively. The results obtained by seven interval robust optimization methods are shown in Table 8. It is observed that (1) for feasibility, all the methods can find feasible solutions and satisfy the constraints, (2) for robustness, the methods using IK-GA-RO, STA-RO, I-STA-RO, and H-STA-RO can obtain robust solutions since their robustness violation values (R) equal to 0, and (3) for optimality, although the result of the SQP-RO method has smaller objective function value, it violates the robustness requirements; within the robust solutions, the result of the H-STA-RO method is superior to others because of its lower f value. Moreover, the results of p r , FE, and T demonstrate the efficiency of the H-STA-RO method. Figure 13 shows the Monte Carlo test results of the deterministic solution and the robust solution of the H-STA-RO method. The deterministic optimum become infeasible in some cases, but the robust optimum of the H-STA-RO method is always feasible even when the decision variables are subject to interval variations.

Example 7: Power Scheduling Design
Power scheduling optimization is an important issue for energy consumption in the process industry. Take the zinc electrowinning process as an example: it accounts for 80% of the total energy consumption of zinc hydrometallurgy. Based on the power time-of-use pricing policy, the aim of power scheduling optimization is to minimize the electricity charge (y) by adjusting the current density (x 1 ), the concentration of Zn 2+ (x 2 ), and H + (x 3 ) in different periods.
A zinc electrowinning process shown in Figure 14 contains seven series potrooms and each potroom has several parallel electrolytic cells. With an appropriate current and zinc acid ratio, zinc ions are deposited on the cathode surface. The electricity use is charged at different prices during three different time periods (peak, shoulder, and off-peak); it is a practice to produce more with lower electricity price. To analyze the power scheduling system, the zinc electrowinning process model is established based on electrochemical reaction mechanism and historical data. With the deterministic parameter, the best result is y = 1.5922 × 10 6 , corresponding to x 1 = [261, 317, 650], x 2 = [60, 45, 60], x 3 = [200, 200, 200]. With incomplete knowledge of the process model, it is more accurate to estimate some parameters as interval values. The uncertain power scheduling optimization problem is formulated as follows:  Figure 14. Zinc electrowinning process.
The parameters in the above problem are as follows: J 0 is the capacity electricity charges; P i is the electricity price in the ith period; T i is the duration of the ith period; V ij is the cell voltage of the jth plant in the ith period; L ij is the magnitude of the current of the electrolysis process of the jth plant in the ith period; g 0 is the zinc daily output; q is the electrochemical equivalent of zinc; E i is the current efficiency in the ith period; N j is the number of cells in the jth plant; B j is the number of plates in a cell in the jth plant; and S is the area of the cathode plate.
The results obtained by five interval robust optimization methods are shown in Table 9. Compared to GA-RO and SQP-RO, the methods based on STA obtain smaller function values corresponding to a more accurate solution. Meanwhile, methods based on STA could obtain a more robust result. Compared to STA-RO and I-STA-RO, H-STA-RO obtains an accurate solution with less function evaluation and runtime, which denotes the efficiency of the proposed method. Figure 15 shows the Monte Carlo test results of the deterministic solution and the robust solution of the H-STA-RO method. The deterministic optimum becomes infeasible in some cases, but the robust optimum of the H-STA-RO method is always feasible.

Conclusions
A hybrid state transition algorithm is proposed to alleviate the problems of robust optimization, including high computational costs and poor convergence. Based on the worst-case analysis, the robust optimization problem can be transformed into a minmax problem. In the outer optimization process (minimization problem), the hybrid state transition algorithm is used to improve the rate of convergence and avoid the local optimum distraction. Meanwhile, the method of sequential quadratic programming is used to strengthen local search ability and reduce computational costs. In the inner optimization process (maximization problem), the second-order Taylor series surrogate model is used to approximate the nonlinear functions and decrease the computational cost. Moreover, to balance the robustness and optimality of candidate solutions, a novel feasibility-checking mechanism is proposed to operate the inner optimization only when a feasible solution is found. Verifying the robustness of the proposed method is conducted using seven examples. The results show that the proposed method offers competitive performance compared with existing robust optimization methods in convergence and efficiency.
In our future work, the robust optimization method for other forms of uncertainties (such as fuzzy uncertainty and interval fuzzy uncertainty) and the applicability enhancement of the surrogate model will be investigated.

Data Availability Statement:
The data presented in this study are available in the verification section of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.