Multi-objective Robust Optimization using a Post-optimality Sensitivity Analysis Technique: Application to a Wind Turbine Design

A general formulation of MOOP is given in Eqn. (1):

where ${\bf f}({\bf x},{\bf p})=[f_{1}~f_{2}~\ldots~f_{m}]^{T}$ denotes the $m$-dimensional vector of performance functions. ${\bf x}=[x_{1}~x_{2}~\ldots~x_{n}]^{T}$ denotes the $n$-dimensional vector of DVs. Note that the nominal values of DV are controllable, $x_{l}^{l}$ and $x_{l}^{u}$ are the lower and upper bounds of $x_{l}$ respectively. ${\bf p}=[p_{1}~p_{2}~\ldots~p_{r}]^{T}$ denotes the $r$-dimensional vector of DEP, which cannot be adjusted by the designer, and they are uncontrollable parameters. The functions $g_{k}({\bf x},{\bf p})=[g_{1}~g_{2}~\ldots~g_{q}]^{T}$ are the constraints. A design that does not violate any constraint is called “feasible”. The set of feasible solutions is called feasible set and named $\mathcal{F}$.

Since there are some trade-offs amongst the $m$ conflicting objectives, none solution of $\mathcal{F}$ can dominate another solution of $\mathcal{F}$ . The optimization problem (1) generally has more than one optimal solutions. Those solutions are defined as Pareto optimal solutions, which cannot be dominated by any other feasible solution [38, 39, 40, 26]. The set of all Pareto optimal solutions is called Pareto optimal set and named $\mathcal{P}$. The Pareto optimal solutions lie on a boundary in the Performance Function Space (PF-Space), called Pareto front.

The concept of robustness is used in many fields such as engineering, biology, economy, computer science [41, 42, 26, 43, 44, 45, 46, 47]. In this paper, the robustness of a product is defined as follows:

Robust design, a term originally introduced by Genichi Taguchi [48], is a way of improving the quality of a product by minimizing the effect of variations, without eliminating the causes themselves. Many researchers refer to robustness [49, 50, 51, 52, 1, 53, 54, 55, 56].

Although different expressions are used, their meanings are similar. In this paper, we define the robust design as follows:

A MOROP aims to find out a solution that is feasible, optimal and robust. To come up with an feasible, optimal and robust design, the three following scenarios have been identified.

1. 1.

   Optimal design and robust design are equally important: It means that performance functions and robustness functions are optimized simultaneously. The designer can use the new performance functions by adding the effects of robustness, instead of original performance functions [3, 4, 6, 7, 57]. Moreover, the robustness functions can be considered as new constraints [8, 9, 58, 59, 14, 16, 17, 60] or additional optimization objectives [15, 61, 62, 63, 64, 65, 66, 67].

2. 2.

   Optimal design is primary and robust design is secondary: It is one post-optimality approach. The final solution is selected from the Pareto optimal set based on the robustness criterion [68, 69, 70, 24, 25, 71, 18, 19, 72].

3. 3.

   Robust design is primary and optimal design is secondary: The final solution is selected from the most robust solutions based on the values of the performance functions [12].

In this paper, the optimal design is supposed to be of primary importance and the robust design of secondary importance. Although there exist some papers dealing with this method  [68, 69, 70, 24, 25, 71, 18, 19, 72], which is also called as post-optimality sensitivity analysis, the final solution is selected from the Pareto optimal set based on its robustness. Moreover, in this paper we consider both small variations in design parameters and large variations in DEPs simultaneously.

Figure 1 illustrates the solutions of a simple multi-objective optimization problem. It has two performance functions ($f_{1}$,$f_{2}$,), some DVs ($\bf x$) and DEPs ($\bf p$). Let us select solution A, solution B and solution C, which are located on the Pareto front. When there are small variations in DVs and DEPs, then their performances varies in the grey area around their nominal values in the PF-Space, as shown in Fig. 1. Obviously, in the PF-Space, solution B has smaller variations around its nominal values than solutions A and C. It means than solution B is more robust than solutions A and C. However, the size of the grey areas associated with solution A and solution C are the same, but their shapes are different. So the comparison between the robustness of solution A and solution C is not an easy task. Toward the different types of uncertainties, there are different methods to define RI for MOOP, such as using the worst case scenario [14, 15, 16, 73, 17], the expectancy measure  [6, 7, 3, 4, 5, 8, 9, 49, 2, 74], and the probabilistic threshold method [8, 9, 2].

We select a RI against the small variations in DVs and DEPs, based on the expectancy measure. A new RI ($I_{RS}$) is introduced and is based on the actual variations of performances. Under the small variations in DVs and DEPs, the actual performances are distributed around the nominal values. In the PF-Space, the actual performances follow a multivariate distribution. Here, in order to simplify the index associated with each performance function, the standard deviation $\sigma$ of the actual performances is used as a measure for the robustness: the smaller the standard deviation, the more robust the design. The absolute value of the difference between the expected value ($\mu$) and the nominal value ($f_{0}$) is also a robustness measure: the smaller the absolute difference, the more robust the design. Figure 2 shows an example, which follows a normal distribution, $\sigma_{f_{i}}$ and $\mu_{f_{i}}$ are the standard deviation and expected value of the $i$th performance function under small uncertainties; $f_{i,0}$ is the nominal value; $f^{max}_{i}$ and $f^{min}_{i}$ are the maximum and the minimum values of the $i$th performance function on the Pareto front, respectively.

Here, we assume that the standard deviation ($\sigma$) and the absolute value of the difference between $\mu$ and $f_{0}$ have the same importance for the designer. Moreover, we assume that the robustness of all performance functions has the same importance. The robustness index $I_{RS}$ of a MOOP with respect to small variations in DVs and DEPs is defined as a scalar. The robustness of each performance function is also normalized. To normalize the sum of the standard deviation ($\sigma_{f_{i}}$) and the absolute difference ($|\mu_{f_{i}}-f_{i,0}|$), we divide it by the difference between the two extreme values of the $i$th performance function, namely, $f^{max}_{i}-f^{min}_{i}$. As a result, $I_{RS}$ is defined as follows:

The smaller $I_{RS}$, the more robust the design.

Note that even if the variations are small, the constraints may not be satisfied due to variations. This refers to the robustness of constraints (reliability). As a matter of fact, the index $I_{RS}$ is not discussed deeply in this paper. It is just based on the distributions of the performance function values with respect to these small variations. To simplify the index, the reliability with regard to small variations is not considered in this paper.

Since the DVs are controllable and the DEPs are uncontrollable, we assume that there are only large variations in DEPs. The mapping functions between the DVs and performance functions can change a lot due to large variations in DEPs. For a MOOP while considering large variations in DEPs, we assume that the initial DEPs are ${\bf p}={\bf p}_{0}$, the set of feasible solutions is named $\mathcal{F}_{0}$ and the set of Pareto optimal solutions is denoted $\mathcal{P}_{0}$. It is noteworthy that those Pareto optimal solutions are alternative solutions for the designer.

Since large variations in DEPs exist, the DEPs may change from ${\bf p}_{0}$ to ${\bf p}_{new}$. The design environment parameters $\bf p$ are supposed to take $N$ discrete values: ${\bf p}_{1},{\bf p}_{2},\cdots,{\bf p}_{N}$. The Probability Density Functions (PDF) of $\bf p$ are $h({\bf p}_{1}),h({\bf p}_{2}),\cdots,h({\bf p}_{N})$. The initial DEPs ${\bf p}_{0}$ are equal to the ones having the maximum PDF amongst the $N$ discrete values. The corresponding feasible sets are $\mathcal{F}_{1},\mathcal{F}_{2},\cdots,\mathcal{F}_{N}$. A feasible and Pareto optimal solution in $\mathcal{P}_{0}$ may not be Pareto optimal, and not feasible in the new environments. To compare the solution’s robustness against large variations in DEPs, the traditional methods are not applicable. Toward a MOOP, a definition of the solution’s robustness against large variations in DEPs is proposed thereafter:

Figure 3 illustrates the solutions of a simple MOOP. Similarly to Fig. 1, the problem has two objective functions ($f_{1}$,$f_{2}$), some DVs ($\bf x$) and DEPs ($\bf p$). ${\bf p}_{0}$ is the nominal values of DEPs. Let us consider Pareto-optimal solutions $A$ and $B$ as examples. The grey area shows the variations in performances of solution $A$ and solution $B$ in the PF-Space, under small variations in DVs and DEPs. As seen before, we can conclude that solution $B$ is more robust than solution $A$ with regard to small variations in DVs and DEPs.

On the contrary, it may not be the case if there are large variations in DEPs. As shown in Fig. 3, in a new environment, we assume that when the DEPs ${\bf p}$ change from ${\bf p}_{0}$ to ${\bf p}_{new}$, there are large variations in the PF-Space for all the solutions. For instance, in the PF-Space, solution A moves from $A_{0}$ to $A_{new}$, solution B moves from $B_{0}$ to $B_{new}$. Both solutions have quite large variations in $\bf f$. Then, the following question remains: “How can we compare the robustness of solutions $A$ and $B$”?.

As shown in Fig. 3, toward solution $A$, $\Delta f_{1,A}^{S}$ and $\Delta f_{2,A}^{S}$ represent the largest distance of the actual performance function values and nominal performance function values in $f_{1}$ and $f_{2}$ respectively, with regard to the small variations in DVs and DEPs. $\Delta f_{1,A}^{L}$ and $\Delta f_{2,A}^{L}$ represent the largest distance of the actual performance function values (assuming only two samples: $A_{0}$ and $A_{new}$) and nominal performance function values in $f_{1}$ and $f_{2}$ respectively, with regard to large variations in DEPs. If $\Delta f_{1,A}^{L}>>\Delta f_{1,A}^{S}$, $\Delta f_{2,A}^{L}>>\Delta f_{2,A}^{S}$ , and $\Delta f_{1,B}^{L}>>\Delta f_{1,B}^{S}$, $\Delta f_{2,B}^{L}>>\Delta f_{2,B}^{S}$, we can see that the caused actual performance function values of both solution $A$ and solution $B$ are sufficiently far from their nominal ones, with regard to large variations in DEPs. As a consequence, we can conclude that the traditional methods, which provide some results based on the difference between actual performance function values and nominal ones, make little sense. Therefore, the method proposed in this paper aims to compare their relative positions in the PF-Space associated with the new environment. From Figure 3, we can see that, in the new environment, ${\bf p}={\bf p}_{new}$, solution $A$ is still on the new Pareto front, which means that it is still one of the best solutions for the designer. Meanwhile, solution $B$ is far away from the new Pareto front, which means that it is no longer a good choice for the designer. As a result, we can conclude that solution $A$ is more robust than solution $B$ with regard to the large variations in DEPs.

As a result, a RI with regard to large variations in DEPs is defined thereafter. Toward the discrete probability distribution of DEPs, mathematically, the RI $I_{RL}$ of a solution ${\bf x}$ is defined as follows:

where $I_{F}$ is the Feasibility Index of the solution; $I_{P}$ is the Pareto optimality Index of the solution; $h(\bf p)$ is the PDF of $\bf p$; $I_{rank}({\bf x},{\bf p}_{j})$ is the individual’s ranking in the new environment where ${\bf p}={\bf p}_{j}$ and amounts to the number of individuals by which it is dominated amongst the alternative solutions, plus one [75, 76]. For a better understanding, a simple example is shown in Fig. 4. In a new environment, if a solution is still non dominated by any other solution, then the $I_{P}$ value will be equal to one for that solution. Otherwise, the $I_{P}$ value will be lower than one, but greater than zero. Note that the number of the alternative solutions affects the value. However, the proposed definition can divide the alternative solutions into different groups based on their robustness.

Note that the index $I_{RL}$ is bounded between zero and one. On the one hand, if a solution is feasible in all environments and cannot be dominated by any other solution in all possible environments, then $I_{RL}=0$. On the other hand, if a solution is non-feasible in some new environments, then $I_{RL}=1$. If each solution belongs to the set of Pareto optimal solutions $\mathcal{P}_{0}$, then if it is feasible in all environments and its $I_{RL}$ value will be greater than or equal to zero and smaller than one.

Note that the individual’s ranking $I_{rank}({\bf x},{\bf p}_{j})$ corresponds to a specific value of DEPs: ${\bf p}_{j}$, where $j=1,2,\cdots,N$. In case there exist continuous probability distributions of the DEPs, $N$ becomes infinite and it is difficult to assess the robustness index $I_{RL}$ for a solution ${\bf x}$. However, since the domain of the DEPs can be partitioned into many small parts, we can simplify such a problem by using a discrete probability distribution of the DEPs.

Finally, the smaller $I_{RL}$, the more robust the design with regard to large variations in DEPs.

In this paper, we consider not only the RI against small variations, $I_{RS}$, but also the RI against large variations, $I_{RL}$. The RI of a Pareto-optimal solution is represented as a vector. The designer can analyze the robustness of Pareto optimal solutions in the Robustness Function Space (RF-Space). For a better understanding, we take $I_{RS}$ as one dimension of the RF-Space, another dimension of the RF-Space is $I_{RL}$, as shown in Fig. 5.

Thanks to the proposed method, each Pareto optimal solution has a corresponding position in the RF-Space, as shown in Fig. 5. If $I_{RS}$ and $I_{RL}$ are not conflicting, then the designer will be able to select the most robust solution immediately, namely, the solution that minimizes both $I_{RS}$ and $I_{RL}$. If $I_{RS}$ and $I_{RL}$ are two conflicting objectives, then a new Pareto front in the RF-Space will appear. The Pareto-robust solutions represent the most robust solutions amongst the Pareto optimal solutions. The set of Pareto-robust solutions is named $\mathcal{P}_{R}(\mathcal{P})$. Finally, the designer can select the final solution from this set according to his/her requirements.

A flow chart is illustrated in Figure 6 for a better understanding of the proposed post-optimality sensitivity analysis technique. The first step is to determine whether the design environment parameters (DEP) are subject to large variations or not. The distinction between small and large variations is made with regard to their effect on the performance functions. The small variations in design parameters refer to the ones that have a small effect on the performance functions. The large variations in design parameters refer to the ones that have large effect on the performance functions. In the second step, Pareto optimal solutions are obtained for different values of DEP. Those different values of DEP are due to large variations in DEP. In the third step, robustness indices $I_{RS}$ and $I_{RL}$ are calculated for each solution. The new Pareto ranking of the obtained Pareto optimal solutions are obtained in different design environments. Finally, the designer can make a decision from the obtained Pareto optimal solutions based on their robustness, namely, based on the values of $I_{RS}$ and $I_{RL}$ for each solution.

The problem is defined in Eqn. (4). There are two performance functions $f_{1}$ and $f_{2}$, one design variable $x$, one design environment parameter $p$ and one constraint:

$x$ is subject to small variations and the variations in $x$ follow a uniform distribution and are bounded between $-0.1$ and $0.1$. $p$ is subject to large variations and can take three possible values: $p=3$, $p=5$ and $p=8$. The probabilities for $p$ to take those three values are $h(1)=0.2$, $h(2)=0.5$ and $h(3)=0.3$, respectively. The initial value of $p$ is equal to 5.

Figure 7 shows the values of the performance functions $f_{1}$ and $f_{2}$ for $p=3$, $p=5$ and $p=8$, respectively, i.e., in different design environments. For $p=5$, let us consider the five Pareto optimal solutions: solution A obtained for $x=1$; solution B obtained for $x=2$; solution C obtained for $x=3$; solution D obtained for $x=4$ and solution E obtained for $x=5$. These solutions form the Pareto front in the PF-Space and are alternative solutions for the designer. Their new positions for $p=3$ and $p=8$ in the PF-Space are shown in Fig. 7. The vertical line represents the constraint related to the first performance function. Note that the solutions to the right of this line are feasible, whereas the solutions to the left are not feasible.

To assess the robustness index $I_{RS}$ for each alternative solution, the samples around the nominal values are generated by Latin Hypercube Sampling (LHS) [77, 78, 79]. 1000 samples for each alternative solution are generated with regard to the small variations in DV. Then $I_{RS}$ for each alternative solution can be assessed with Eqn. (2).

To determine the robustness index $I_{RL}$ associated with each alternative solution, their new positions should be located when $p$ changes from its initial value to its new value. The results can be found from Fig. 7. When $p=3$, solution A becomes unfeasible, then the index $I_{F}$ defined in Eqn. (3d) is null for solution A. Solution B and solution C become Pareto optimal amongst these five solutions. Solution D and Solution E are dominated by other alternative solutions. When $p=8$, all the alternative solutions are feasible and can not be dominated by others. Then $I_{RL}$ for each alternative solution can be calculated with Eqn. (3).

Note that all Pareto optimal solutions have sufficiently large variations in performance functions with regard to large variations in DEPs. In Tab. 2, $f_{1}$ and $f_{2}$ denote the nominal values of the performance functions for the alternative solutions. $\Delta f_{1}^{S}$ and $\Delta f_{2}^{S}$ represent the largest distance of the actual performance function values (1000 samples for each alternative solution) and nominal performance function values in $f_{1}$ and $f_{2}$ respectively, with regard to the small variations in DV. $\Delta f_{1}^{L}$ and $\Delta f_{2}^{L}$ represent the largest distance of the actual performance function values (3 samples for each alternative solution) and nominal performance function values in $f_{1}$ and $f_{2}$ respectively, with regard to large variations in DEP. If the traditional methods (results based on the difference between actual performance function values and nominal ones) are selected to measure the robustness of the alternative solutions, here is the order of the six solutions from the most robust one to the least robust one: (1) E; (2) D; (3) C; (4) B; (5) A. Even if solution E is the most robust one, $\Delta f_{1,E}^{L}$ and $\Delta f_{2,E}^{L}$ are still quite large. Obviously, $\Delta f_{1,E}^{L}>>\Delta f_{1,E}^{S}$ and $\Delta f_{2,E}^{L}>>\Delta f_{2,E}^{S}$, so with regard to large variations in DEPs, the traditional methods make little sense.

Here is the order of the six solutions from the most robust one to the least robust one while using the robustness index $I_{RL}$ defined in Eq. (3a): (1) B; (1) C; (3) D; (4) E; (5) A. Solution B and solution C are Pareto optimal no matter the value of DEP $p$. Therefore, the designer may select solution B or solution C instead of solution E.

Figure 8 illustrates the positions of the alternative solutions in the RF-Space. In this space, one dimension is $I_{RS}$, another dimension is $I_{RL}$. In this example, $I_{RS}$ and $I_{RL}$ are two conflicting objectives, as shown in Fig. 8. Solutions C, D, E form a new Pareto front in the RF-Space, which represents the most robust solutions amongst these five alternative solutions. For a better understanding, the performance function values and their RIs are given in Tab. 2. The indices $I_{F}$ and $I_{P}$ are also given in Tab. 2.

Thanks to the proposed method, the designer can select the final solution from the new Pareto front in the RF-Space, according to his/her requirement. For example, if the designer prefers the $I_{RS}$ index, solution E will be selected. If the designer prefers the $I_{RL}$ index, solution C will be selected.

According to the proposed method, the final solution comes from the Pareto optimal set. In this paper, the Pareto optimal solutions are obtained by using the genetic algorithm (NSGA II  [81]). The obtained Pareto front $\mathcal{P}$ for this problem is illustrated in Fig. 10, including 200 alternative solutions.

To assess the robustness index $I_{RS}$ for each alternative solution, the samples around the nominal values are generated by Latin Hypercube Sampling (LHS) [77, 78, 79]. 1000 samples for each alternative solution are generated with regard to the small variations in DVs and DEPs. Then $I_{RS}$ for each alternative solution can be assessed using with Eqn. (2).

To determine the robustness index $I_{RL}$ with regard to large variations in one DEP for each alternative solution, the new positions of the alternative solutions should be presented, when DEPs change from initial values to the new values. Then the feasibility index and the Pareto optimality Index of a solution can be determined.

In this problem, the wind speed has a continuous probability distribution, varying from 6 m/s to 14 m/s. To simplify the problem, we calculate the probability distribution of the wind speed, and convert it into a discrete probability distribution problem. Table 5 shows the probability distribution of the wind speed. The solution having the maximum PDF amongst the N probabilities is appointed as the initial DEPs ${\bf p}_{0}$, hence ${\bf p}_{0}={\bf p}_{5}$ in this problem.

Figure 11 illustrates the positions of the alternative solutions in different design environments. The vertical lines represent the constraints in the first objective function, i.e., the produced power $P$, then the solutions between the two lines are feasible. From the results, we can see that some alternative solutions are non-feasible in the new environments. They are illustrated by red crosses in Fig. 11. The index $I_{F}$ defined in Eqn. (3d) is null for those solutions. In the new environment, some alternative solutions are dominated by other alternative solutions, then their index $I_{P}$ defined in Eqn. (3e) is lower than one, but greater than zero. Then the $I_{RL}$ for each alternative solution can be calculated with Eqn. (3a).

Figure 12 illustrates the positions of the alternative solutions in the RF-Space. In this space, one dimension is $I_{RS}$, another dimension is $I_{RL}$. In this example, $I_{RS}$ and $I_{RL}$ are two conflicting objectives, as shown in Fig. 12. A new Pareto front in the RF-Space appears, which represents the most robust solutions amongst the Pareto optimal solutions and denoted by $\mathcal{P}_{R}(\mathcal{P})$.

For a good understanding of the two proposed robustness indices defined by Eqn. (2) and Eqn. (3a), the six solutions A, B, C, D, E and F are selected. Their corresponding positions in the PF-Space and the RF-Space are illustrated in Fig. 10 and Fig. 12. All alternative solutions and the initial design are shown in Fig. 13. The initial design, solution F, is plotted in black and the other five solutions are plotted in yellow.

Solution A has the minimum $I_{RS}$ value in the $\mathcal{P}_{R}(\mathcal{P})$. Solution B dominates Solution F in the RF-Space. Solution C has the minimum $I_{RL}$ value in the $\mathcal{P}$. Solution D can be dominated by many solutions such as Solutions B, C and E in the RF-Space. Solution E dominates the initial design in the PF-Space. The 3D models of the six selected solutions and the values of their DVs are depicted in Fig. 14.

Note that the value of index $I_{P}$ of a solution depends on its ranking amongst the Pareto optimal solutions. To determine the index $I_{P}$ of the initial design, i.e., solution F, we include it into the Pareto optimal set (including 200 alternative solutions) and rank the 201 individuals to generate its ranking in different DEPs. The performance function values and the RIs for these six solutions are given in Tab. 6. For a better understanding of index $I_{RL}$, indices $I_{F}$ and $I_{P}$ of these six solutions are also given in Tab. 6.

With regard to the small variations in DVs and DEPs, each solution has 1000 generated samples in the PF-Space. Figure 15 shows the samples and the nominal values of these six solutions. It is apparent that Solution A is the most robust one amongst the six solutions considering small variations in DVs and DEPs as the dispersion of the two performance functions is a minimum with regard to small variations in DVs and DEPs for this solution. It matches with their positions in the RF-Space. The robustness indices $I_{RS}$ with regard to small variations for the five other solutions are given in the fourth column of Tab. 6. Here is the order of the six solutions from the most robust one to the least robust one with respect to $I_{RS}$: (1) A; (2) B; (3) F; (4) E; (5) C; (6) D.

With regard to large variations in DEPs, all Pareto optimal solutions have sufficiently large variations in performance functions. In Tab. 7, $-P$ and $F_{a}$ represent the nominal performance function values of the alternative solutions. $\Delta P^{S}$ and $\Delta F_{a}^{S}$ represent the largest distance of the actual performance function values (1000 samples for each alternative solution) and nominal performance function values in $P$ and $F_{a}$ respectively, with regard to the small variations in DVs and DEPs. $\Delta P^{L}$ and $\Delta F_{a}^{L}$ represent the largest distance of the actual performance function values (9 samples for each alternative solution) and nominal performance function values in $P$ and $F_{a}$ respectively, with regard to large variations in DEPs. The results show that $\Delta P^{L}>>\Delta P^{S}$ and $\Delta F_{a}^{L}>>\Delta F_{a}^{S}$, so with regard to large variations in DEPs, the traditional methods (results based on the difference of actual performance function values and nominal ones) make little sense. On the contrary, the method proposed in this paper is relevant. Indeed, here is the order of the six solutions from the most robust one to the least robust one based on the $I_{RL}$ index: (1) C; (1) B; (3) E; (4) D; (5) F; (6) A.

In summary, the designer can select the final solution from the new Pareto front in the RF-Space, according to his/her requirement. For example, if the designer prefers the $I_{RS}$, solution A should be selected. If the designer prefers the $I_{RL}$, solution C should be selected.