A Crossrate-Based Approach for Reliability-Based Multidisciplinary Dynamic System Design Optimization

Abstract

In practical applications, the multidisciplinary dynamic system design optimization (MDSDO)-based solution is limited by uncertainty, which causes random variation in the physical design variable in the static discipline and the equation of state in the dynamic discipline. To address the lack of reliability of the MDSDO solution, a crossrate-based MDSDO approach (C-MDSDO), consisting of the MDSDO stage and a reliability assessment stage, is proposed in this paper. In the reliability assessment stage, a sub-optimization problem based on the crossrate of the objective reliability index sample trajectory is designed to obtain the shifting vector, which is employed to obtain a sufficiently reliable solution. In addition, the proposed approach adopts a sequential problem-solving framework that avoids nested optimization and a reliability assessment. One numerical case and two engineering cases were employed to validate the effectiveness of the proposed method. The results show that the reliability of the proposed solutions significantly improved.

1. Introduction

In engineering applications, designing a complex system such as a flying vehicle always involves time-independent static and time-dependent dynamic disciplines. In the traditional design method, the design process of the physical design of static disciplines and the control strategies of dynamic disciplines are performed separately and sequentially. The control strategy is developed based on the physical design solution, which is completed beforehand [1,2]. In the design of static disciplines, the design variables involved in different disciplines are highly coupled, which makes traditional design approaches inefficient. Therefore, researchers proposed the multidisciplinary design optimization (MDO) theory, which decouples the various disciplines and significantly improves problem-solving efficiency [3]. The commonly used MDO approaches can be categorized into two types: the single-level approach and the multilevel approach [4]. Commonly used single-level optimization approaches can be further divided into multidisciplinary feasible (MDF) and single-disciplinary feasible (IDF) approaches [5,6]. In contrast, commonly employed multilevel optimization approaches include collaborative optimization (CO), concurrent subspace optimization (CSSO), bi-level integrated system collaborative optimization (BLISS), and analytical target cascading (ATC) [7,8,9,10]. In order to improve the applicability of the above MDO method, researchers propose many improvements [11,12]. As for the design approach for dynamic disciplines, researchers focus on the trajectory optimization of control objects, which is known as the optimal control problem (OCP) [13]. The most frequently used OCP approaches can be divided into indirect and direct approaches [14]. In the indirect approach, the boundary condition and the Hamilton–Jacobi information of the control object are derived on the basis of the maximum principle, and the control problem is converted into a boundary value problem [15,16]. In the direct approach, i.e., direct transcription (DT), the control problem is converted into nonlinear programming (NLP) using discrete control and state variables at the time nodes. Typical direct approaches include control vector parameterization and orthogonal collocation [17,18].

The design approaches for the static and dynamic disciplines are fairly well developed and enable the optimal solution to be obtained for each discipline. However, it is difficult to obtain the optimal solution for a whole system with a split design process [9]. Therefore, researchers proposed the multidisciplinary dynamic system design optimization (MDSDO) theory by introducing the concept of integrated design. Azad et al. coupled the optimization process of static and dynamic disciplines to form the nested optimization algorithm [19]. Furthermore, various researchers constructed the surrogate model of the right-hand-side function of the state equation in the inner dynamic discipline optimization loop, which significantly improves the efficiency of the inner loop and overcomes the inefficiency of the nested algorithm [20,21,22]. Moreover, other researchers employed the simultaneous optimization approach [23,24]. The research on MDSDO is still in its developmental stage, and its problem-solving approaches need to be further improved.

In the static discipline of an MDSDO problem, design variables may vary because of certain known or unknown factors, which is known as the “uncertainty”. To avoid the negative effects associated with uncertainty in design approaches using static disciplines, researchers proposed reliability-based design optimization (RBDO) [25,26], which assesses the deterministic solution and improves its reliability based on the assessment result. In addition, the practicality of the approach is significantly enhanced by further improvements [27,28]. Since the static discipline is an essential part of MDSDO, introducing the reliability assessment also acts to improve reliability-based design optimization for the MDSDO problem. The findings related to reliability-based MDSDO remain limited. Some researchers propose a three-loop framework, in which the reliability assessment is directly nested in the static discipline design loop [29]. However, as the complexity of the problem increases, the efficiency of this nested framework becomes insufficient. Certain researchers propose a transformation-based reliability-based MDSDO framework by transforming the uncertainty design variables into equivalent deterministic and concomitant variables and adding equality constraints in the optimization process of the dynamic discipline [30]. However, this approach changes the original control problem of the dynamic discipline, and the optimization results in the literature show that the reliability indexes of the obtained solutions do not reach the necessary reliability requirements.

In this paper, a crossrate-based MDSDO approach (C-MDSDO) is proposed, which avoids the nesting of the reliability assessment loop and can stably obtain a solution that meets the reliability requirements. C-MDSDO is a single-loop framework that consists of a deterministic MDSDO stage and a reliability assessment stage. In each iteration, the obtained solution from the deterministic MDSDO stage is evaluated in the reliability assessment stage, and the shifting vector that guides the deterministic solution to meet the reliability requirement is obtained according to the assessment result. Since the shifting vector in reliability-based MDSDO must simultaneously consider the impact of uncertainty on the design variables, the state variables, and the control variables, the traditional criteria for the vector solution in the traditional SORA is not applicable. Therefore, a novel sub-optimization problem is formulated. In the formulated problem, the uncertainty samples of the design variables that satisfy the objective reliability index are employed as the optimization object to obtain the inverse most probable point (iMPP)-based shifting vector, which is employed to obtain the initial values of the physical design variables for the next iteration. As a result of the proposed formulation, the difficulty associated with searching the shifting vector due to the decoupling of state variables, control variables, and design variables in the constraints is resolved. Furthermore, our method increases the vector searching efficiency when multiple MPPs exist in highly nonlinear constraints. In the proposed approach, the properties of the original dynamic discipline are maintained, and the framework structure is simple. Several numerical and engineering cases were used to verify the effectiveness of the proposed approach. Moreover, the results from various studies were used for comparison. The results demonstrate the effectiveness and reliability of the proposed method.

The rest of this paper is organized as follows. In Section 2, the definition of the reliability-based MDSDO problem in this paper is given. In Section 3, details of the construction and the process used for the sub-optimization problem are described. In addition, the algorithm flow of the proposed C-MDSDO is elaborated. In Section 4, several case studies are presented that verify the effectiveness and reliability of the proposed approach. The test results are discussed in Section 5, and conclusions are drawn.

2. Multidisciplinary Dynamic System Design Optimization and Uncertainty

In the design of a complex system, the coupling of static and dynamic disciplines is almost unavoidable. MDSDO theory achieves the unification of two types of discipline design optimization through decomposition and coordination. An MDSDO problem can be expressed as the following mathematical expression:

$$\begin{array}{cc}\hfill \underset{d,x\left(t\right),u\left(t\right)}{min}& \varphi \left(d,x\left(t\right),u\left(t\right),t\right)\hfill \\ \hfill subjectto& \dot{x}-f\left(d,x\left(t\right),u\left(t\right),t\right)=0\hfill \\ & h\left(d,x\left(t\right),u\left(t\right),t\right)=0\hfill \\ & g\left(d,x\left(t\right),u\left(t\right),t\right)\ge 0\hfill \end{array}$$

(1)

where $\varphi (·)$ is the optimization objective, $\dot{x}(·)$ is the state equation, $h(·)$ is the equality constraint, $g(·)$ is the inequality constraint, $d$ represents the design variable, $x\left(t\right)$ represents the state variable, and $u\left(t\right)$ represents the control variable.

An MDSDO problem is commonly solved using the co-design framework in MDF, which is a two-loop structure flow in which the outer loop is design variable optimization and the inner loop is control optimization [31]. The flow chart of the co-design framework is shown in Figure 1.

The outer loop of static discipline design optimization can be expressed as follows:

$$\begin{array}{cc}\hfill \underset{d}{min}& {\varphi }^{*}\left(d,x\left(t\right),u\left(t\right),t\right)\hfill \\ \hfill subjectto& g\left(d,x\left(t\right),u\left(t\right),t\right)\ge 0\hfill \end{array}$$

(2)

where the state variable $x\left(t\right)$ and the control variable $u\left(t\right)$ are derived from the inner loop.

The optimization problem of the inner loop is expressed in the following, where the design variable $d$ is provided by the outer loop.

$$\begin{array}{cc}\hfill \underset{x\left(t\right),u\left(t\right)}{min}& \varphi \left(d,x\left(t\right),u\left(t\right),t\right)\hfill \\ \hfill subject to & \dot{x}-f\left(d,x\left(t\right),u\left(t\right),t\right)=0\hfill \\ & h\left(d,x\left(t\right),u\left(t\right),t\right)=0\hfill \\ & g\left(d,x\left(t\right),u\left(t\right),t\right)\ge 0\hfill \end{array}$$

(3)

As mentioned in the introduction, the internal optimization problem expressed above is also known as the OCP. The OCP is commonly solved by DT, which transforms the infinite-dimensional original dynamic problem into a finite-dimensional NLP using discrete variables on the time interval. Sequential quadratic programming (SQP) is employed to solve the NLP as expressed below.

$$\begin{array}{cc}\hfill \underset{{{\rm M}}_x,M_u}{min}& \varphi \left(d,{{\rm M}}_x,M_u\right)\hfill \\ \hfill subject to & \xi \left(d,M_x,M_u\right)=0\hfill \\ & h\left(d,M_x,M_u\right)=0\hfill \\ & g\left(d,M_x,M_u\right)\ge 0\hfill \end{array}$$

(4)

where ${{\rm M}}_x$ represents the discrete state variable matrix, $M_u$ represents the discrete control variable matrix, and $\xi$ represents the discrete equation of state vectors.

Using the co-design framework, the problem of the traditional split sequential approach not being able to obtain a system-level optimal solution was solved by researchers. However, as a result of uncertainty in manufacturing, certain physical design parameters cannot reach an ideal value. Under such conditions, the performance of static disciplines is affected, and the state trajectory of dynamic disciplines coupled with design variables changes drastically. The continued application of the deterministic MDSDO solution may lead to unexpected failures. Therefore, it is necessary to introduce reliability assessment to evaluate a deterministic MDSDO solution and obtain a reliability-based MDSDO solution.

In this paper, a test case was employed to demonstrate the impact on the dynamic discipline state trajectory when the design variables in the static discipline are subject to uncertainty. The Monte Carlo sampling (MCS) technique was used to generate 10,000 MCS samples that obeyed the random distribution of the design variables. These random samples represented the value of the physical variable that could be taken when subject to uncertainty. With the assistance of the known differential equation, the state trajectory of each random sample was inverted. The trajectory of the deterministic MDSDO solution and the trajectories of the corresponding uncertain samples are shown in Figure 2. First, we randomly selected uncertainty samples with deviations of 1σ, 2σ, and 3σ, and their corresponding state trajectories are shown in the figure. As can be observed from the figure, the difference between the corresponding state trajectory and the deterministic state trajectory gradually increases with an increase in deviation of the uncertainty sample. Then, by employing the MCS technique to generate a large number of uncertainty samples, the range of possible trajectories of the deterministic MDSDO design under the impact of the uncertainty can be obtained. This variation range is shown by the green trajectory curves.

The constraint function in this case required the state variable to be no less than −0.4. In Figure 2, there are a large number of green lines crossing the critical value of −0.4, which means that there is a considerable probability of failure when the MDSDO solution is affected by uncertainty. In practical applications, this MDSDO design is clearly not feasible.

This section introduces the definition of MDSDO and discusses the impact of uncertainty on the MDSDO solution, which illustrates the necessity of introducing reliability assessment. In the next section, the definition of the reliability-based MDSDO problem solved in this paper is given, and the proposed C-MDSDO approach is described.

3. The Proposed C-MDSDO

There is no uniform definition of the reliability-based MDSDO problem. Hence, this paper attempts to define the problem and propose the approach from the perspective of practical applications. Under actual circumstances, uncertainty is inherent in both static and dynamic disciplines. However, for the dynamic discipline with equality constraints, considering the constraint to be constant makes the problem complicated or even unsolvable in any case. As an attempt to introduce reliability into MDSDO, the control strategy of the dynamic discipline is regarded as stable and reliable in this paper. At the same time, the design variables of the static discipline suffer from uncertainty. In general, this paper aims to obtain a generalized control strategy that satisfies the constraints of the system when the uncertainty design variable exists.

The reliability-based MDSDO problem solved in this paper is expressed as the following formula:

$$\begin{array}{cc}\hfill \underset{u_d,x\left(t\right),u\left(t\right)}{min}& \varphi \left(u_d,x\left(t\right),u\left(t\right),t\right)\hfill \\ \hfill subject to & \dot{x}-f\left(u_d,x\left(t\right),u\left(t\right),t\right)=0\hfill \\ & h\left(u_d,x\left(t\right),u\left(t\right),t\right)=0\hfill \\ & Prob\left[g\left(u_d,x\left(t\right),u\left(t\right),t\right)\ge 0\right]\ge \Phi \left(\beta \right)\hfill \end{array}$$

(5)

where $u_d$ represents the mean of the design variable, and $\beta$ represents the objective of the reliability index.

As shown in the formula, the design variables here are coupled in dynamic disciplines. As discussed previously, the problem can be solved by discretizing the problem using DT and performing the reliability assessment at each time point. However, the disadvantage in efficiency is clear. Therefore, a sequential C-MDSDO approach is proposed to solve the problem related to the efficiently of the above formula. C-MDSDO consists of a deterministic MDSDO stage and a reliability assessment stage. A sub-optimization problem is constructed to achieve the reliability assessment. In the following section, the overall algorithm flow is presented, and then the constructed sub-optimization problem is described in detail.

3.1. Algorithm Flow of the Proposed Approach

In C-MDSDO, the available co-design algorithm is first employed to obtain a deterministic solution, which is applied for the reliability assessment in the following stage. The reliability assessment is based on a constructed sub-optimization problem, in which the stochastic optimization algorithm is used to obtain the shifting vector. With the assistance of the shifting vector, the constraints of the deterministic MDSDO problem are gradually modified, and the constraints are used for the next iteration. When the objective and reliability index converge, the solution to the MDSDO problem satisfying the reliability index is obtained. This decoupled sequence framework effectively avoids nesting reliability assessments in complex dynamic discipline optimization. The flowchart of C-MDSDO is shown in Figure 3.

The details of each C-MDSDO step are as follows:

(a)

The boundary constraints of the design variables, state variables, and control variables are established given the initial guesses;

(b)

The co-design framework is employed to obtain the solution of deterministic MDSDO of the step $k=1$;

(c)

MCS population is established for calculating the failure probability based on the uncertainty information of design variables. The state trajectories of the MCS population are deduced based on the Runge–Kutta method and the control trajectories $u_k\left(t\right)$. The failure probability is calculated for the state trajectory of the MCS population, and the objective performance is calculated.

In this study, the MCS population number was set to 1e5, and the state trajectory deduction method is described in Section 3.2. Figure 2 in the previous section shows the deduction trajectory $X_M\left(t\right)$ of a test case. In addition, the failure probability and reliability index are calculated as shown in Formula (6).

$$\{\begin{array}{l}\beta ={\Phi }^{-1}\left(pf\right)\\ pf=\frac{num\left[g\left(d_M,X_M\left(t\right),u\left(t\right),t\right)\ge 0\right]}{N_M}\end{array}$$

(6)

where $num\left[g\left(d_M,X_M\left(t\right),u\left(t\right),t\right)\ge 0\right]$ represents the number of MCS populations that violate the constraint.

(d)

If the objective performance is convergence and the reliability index is satisfied, go to g; otherwise, go to e;

(e)

The design variables $d_{\mathrm{k}}$ and control trajectories $u_{\mathrm{k}}\left(t\right)$ in step $\mathrm{k}$ are substituted into the formulated sub-optimization problem for reliability assessment, and the shifting vector $v_k$ is obtained;

(f)

The shifting vector is used to obtain the initial values of the design variables of the next iteration $\mathrm{k}+1$ with $d_{\mathrm{k}+1}=d_k+v_k$. Go back to b;

(g)

The algorithm is terminated, and the reliability-based MDSDO result is output.

This section describes the main framework of the proposed approach and elucidates the relationship between the main problem and the sub-optimization problem. In the next section, the formulation of the sub-optimization problem is described. On the basis of the formulated sub-optimization problem, the problem-solving method of the shifting vector is presented, which is employed to improve MDSDO reliability.

3.2. Reliability Assessment with the Formulated Sub-Optimization Problem

3.2.1. Accompanying Sample for the Formulated Sub-Optimization Problem

The SORA is a popular framework in the static discipline that can solve various reliability-based optimization problems. MPP searching is the core part of the framework, and it determines the calculation of the shifting vector, which is employed to help to find a sufficiently reliable solution. In each iteration step, constraints are shifted with the vector, and the feasible domain gradually shrinks until the reliability index meets the objective requirement.

However, in an MDSDO problem (5), the constraint is usually coupled with the control variable and the state variable. Although it is possible to obtain the i-MPP using the SORAs at each time node through time discretization and calculate the shifting vector using a certain criterion, this will undoubtedly consume a great deal of computational resources. Therefore, in this paper, instead of searching for the i-MPP of the design variable solution, random numbers of the accompanying sample $R$ satisfying the objective reliability index are generated, as shown in Figure 4, where the accompanying samples with the reliability indexes $\beta =1$, $\beta =2$, and $\beta =3$ are shown with different colors.

The generation of the accompanying sample is expressed as in Formula (7):

$$\begin{array}{cc}\hfill find& R\hfill \\ \hfill subjectto& \Vert U_R\Vert ={\beta }_{obj}\hfill \end{array}$$

(7)

where $U_R$ represents the accompanying sample in the normal space, and ${\beta }_{obj}$ is the objective reliability index.

Since the accompanying samples are employed in the sub-optimization problem, excess accompanying samples will increase the burden of the optimization process. Therefore, according to the dimension $dim$ of the reliability-based MDSDO problem, $Nu{m}_R=4dim$ accompanying samples are employed as follows.

(a)

According to the dimension $dim$ of the design variable in the reliability-based MDSDO problem, vectors for the accompanying sample generation are generated. In Formula (8), the values of vector components are all random numbers, and the interval is between 0 and 1;

$${{\rm T}}_R=\left[{\tau }_1,{\tau }_2,{\tau }_3,\dots ,{\tau }_{Nu{m}_R}\right]=\left[\begin{array}{l}{v}_{1,1},v_{1,2},v_{1,3},\dots ,v_{1,dim}\\ v_{2,1},v_{2,2},v_{2,3},\dots ,v_{2,dim}\\ v_{3,1},v_{3,2},v_{3,3},\dots ,v_{3,dim}\\ \dots \\ v_{Nu{m}_R,1},v_{Nu{m}_R,2},v_{Nu{m}_R,3},\dots ,v_{Nu{m}_R,dim}\end{array}\right]$$

(8)

(b)

According to the objective reliability index ${\beta }_{obj}$ of the reliability-based MDSDO problem, the accompanying samples are obtained.

$$R=\frac{{{\rm T}}_R}{\Vert {{\rm T}}_R\Vert }{\beta }_{obj}$$

(9)

3.2.2. Formulated Sub-Optimization Problem and Problem-Solving Process

To avoid neglecting the influence of the objective function performance when optimizing the reliability index of the solution, the objective function performance is introduced into the shifting vector search process to form the suggested sub-optimization problem, as in Formula (10).

$$\begin{array}{ll}\hfill \underset{v^k}{min}& \varphi \left(d_k-v_k,x_{dv}\left(t\right),u_k\left(t\right),t\right)·C_{rate}\left(R_k\right)\hfill \\ \hfill subjectto& \dot{x}-f\left(R_k-v_k,x_{Rv}\left(t\right),u_k\left(t\right),t\right)=0\hfill \\ & d_{min}\le d_k-v_k\le d_{max}\hfill \end{array}$$

(10)

where $k$ represents the iteration step of the MDSDO stage, $v_k$ represents the shifting vector, $d_k$ represents the design variable solution from the MDSDO stage, $u_k\left(t\right)$ represents the control strategy from the MDSDO stage, and $x_{dv}\left(t\right)$ and $x_{Rv}\left(t\right)$ are the deduction trajectories obtained based on $\left(d_k-v_k\right)$ and $\left(R_k-v_k\right)$, respectively. The deduction method is described in the following. $C_{rate}$ is the index of the constraint violation and is denoted the “crossrate”.

$C_{rate}$ counts the number of constraint violations of the trajectory in the observation point, which is deduced based on the accompanying samples $R_k$. Therefore, to illustrate the crossrate calculation method, the deduction process needs to be described.

In the traditional approach, the solution of the design variable $d_k$ from the MDSDO stage is directly used in the reliability assessment stage. The solution of the state and control variable is processed under an ideal situation. However, in a practical situation, the accompanying sample $R_k$ represents the product with the limit reliability index, so that it is on the verge of failure. Without additional calibration, the product $R_k$ will apply the same control strategy $u_k\left(t\right)$ as the ideal product. Therefore, in this paper, the control strategy $u_k\left(t\right)$ that was designed on the basis of the deterministic situation is employed to deduce the state trajectory.

Since the control strategy from the previous process and the equation of state are known, the fourth-order Runge–Kutta method can be adopted for the deduction process of the accompanying samples. This method applies Formula (11) to obtain the scatter of the state trajectory within the given time interval based on the known equation of state, the initial state, and the given step size.

$$\{\begin{array}{l}x\left(m+1\right)=x\left(m\right)+\frac{h}{6}\left(q_1+2q_2+2q_3+q_4\right)\\ x\left(0\right)=x_0\end{array}$$

(11)

where $x_0$ is the initial state, $x\left(m\right)$ are the discrete points at the time node $m$, and $h$ is the given step size. Other parameters are calculated as shown in Formula (12).

$$\{\begin{array}{l}{q}_1=f\left(t\left(m\right),x\left(m\right)\right)\\ q_2=f\left(t\left(m\right)+\frac{h}{2},x\left(m\right)+\frac{h}{2}{q}_1\right)\\ q_3=f\left(t\left(m\right)+\frac{h}{2},x\left(m\right)+\frac{h}{2}{q}_2\right)\\ q_4=f\left(t\left(m\right)+h,x\left(m\right)+h{q}_3\right)\end{array}$$

(12)

where $f(·)$ represents the equation of state of the system and $t(\cdot )$ is the given time node.

At the beginning of trajectory deduction, each accompanying sample is substituted into the original equation of state to obtain the corresponding initial state value. Then, the deduction trajectory of each accompanying sample is obtained with the control strategy $u_k\left(t\right)$ and the given discrete-time $t(·)$ that are obtained in deterministic MDSDO. Once the deduction trajectories of the accompanying samples are obtained, the calculation of the crossrate $C_{rate}$ can be performed as follows.

$$C_{rate}\left(R_k\right)=1+\frac{{Num}_R{}^{\prime }}{{Num}_R}$$

(13)

where ${num}_R$ represents the number of accompanying samples and ${Num}_R{}^{\prime }$ represents the constraint violation number of the accompanying samples.

Thus, the composition of the sub-optimization problem, the deduction of the state trajectory, the selection of the control trajectory, and the calculation of the crossrate are presented. The complete algorithmic flow of the sub-optimization problem is described in the following, and the flow chart of the proposed reliability assessment stage that contains the sub-optimization problem-solving loop is shown in Figure 5.

(a)

In each step $k$, the solution $\left\{d_k,u_k\left(t\right)\right\}$ from the deterministic MDSDO is employed as the initial condition of the sub-optimization problem;

(b)

According to the uncertainty information, the accompanying samples $R_k$ that satisfy the objective reliability index ${\beta }_{obj}$. are generated;

(c)

The random optimization algorithm is employed to solve the sub-optimization problem shown in Formula (10). Particle swarm optimization is employed in this paper;

(d)

In each iteration, according to the control trajectory $u_k\left(t\right)$, the deduction trajectory of the $R_k-v_k$ is obtained. The $C_{rate}$ and the performance of the objective function in Formula (10) of each particle are calculated. The particle optimum and global optimum are compared and updated, the particle position is updated, and the next iteration is performed;

(e)

The convergence condition of the sub-optimization problem is satisfied. The shifting vector $v_k$ is output. Otherwise, the step-d iteration is continued;

(f)

The sub-optimization loop is terminated, and the design $d_k$ in step $k$ is shifted with $d_{k+1}=d_k+v_k$, which is treated as the initial condition of the next deterministic MDSDO.

This section attempts to define the reliability-based MDSDO problem studied in this paper. For definition, the main structure of the C-MDSDO approach proposed in this paper is first introduced, which consists of a deterministic MDSDO stage and a sub-optimization-based reliability assessment stage. Then, the construction of the suggested sub-optimization stage and the algorithm flow are presented. Moreover, the generation of the accompanying samples, the deduction of the state trajectories, and the selection of the control strategies are described. In the next section, several case studies are employed to verify the effectiveness of the proposed C-MDSDO, and the result from a similar study is introduced for comparison.

4. Test Example

Since there is little research on reliability-based MDSDO, it is difficult to use existing results for comparison. In this paper, only the optimization results of the developed algorithm can be utilized in the first case to demonstrate the advantages of the proposed C-MDSDO. Moreover, two other cases were employed to further validate the effectiveness of the proposed approach, one of which is a commonly employed OCP problem in which the physical design variables are usually taken to be constant. The other case is an engineering case of a three-degree-of-freedom robot arm with black-box equations of state.

4.1. Van Der Pol Oscillator

The Van Der Pol Oscillator [30] is a classical numerical case of optimal control with two design variables$\left[{\mathrm{d}}_1{,\mathrm{d}}_2\right]$, two state variables$\left[x_1{,\mathrm{x}}_2\right]$, and one control variable $\left[\mathrm{u}\right]$. In OCP, the design variables are usually taken to be constant values before further trajectory optimization. In this paper, the design variables were also treated as optimization objects and were attributed with uncertainty in order to verify the effectiveness of the proposed approach. The mathematical expression of this case is shown in Formula (14). The original constraint was transformed into a reliability constraint, and the objective reliability index was set as $\beta =3$.

$$\begin{array}{cc}\hfill \underset{d,x\left(t\right),u\left(t\right)}{min}& \varphi ={\displaystyle \underset{0}{\overset{5}{\int }}\left[x_1^2\left(t\right)+x_2^2\left(t\right)+u^2\left(t\right)\right]}dt\hfill \\ \hfill subjectto& {\dot{x}}_1\left(t\right)=x_1\left(t\right)\left(1-x_2^2\left(t\right)\right)-d_1{d}_2{x}_2\left(t\right)+u\left(t\right)\hfill \\ & {\dot{x}}_2\left(t\right)=x_1\left(t\right)\hfill \\ & prob\left[\mathrm{g}=-x_1\left(t\right)-0.4\ge 0\right]\le \Phi \left(-\beta \right)\hfill \\ & \beta =3,{\mathrm{d}}_1\sim N\left({\mu }_{{\mathrm{d}}_1},0.02\right){,\mathrm{d}}_2\sim N\left({\mu }_{{\mathrm{d}}_2},0.03\right)\hfill \\ & \mathrm{d}=\left[{\mathrm{d}}_1{,\mathrm{d}}_2\right],\mathrm{x}=\left[x_1,x_2\right],\mathrm{t}=\left[0,5\right]\hfill \\ & {\mathrm{d}}_{min}\le d\le d_{max}\hfill \\ & {\mathrm{x}}_{min}\le x\le x_{max}\hfill \\ & {\mathrm{u}}_{min}\le \mathrm{u}\le u_{max}\hfill \end{array}$$

(14)

The boundary condition of all the variables is listed in

Table 1. The information of the design, state, and control variable of test 4.1.

Variable	Range		Initial	Final
	Max	Min
$\left[{\mathrm{d}}_1{,\mathrm{d}}_2\right]$	[5, 5]	[−5, −5]	/	/
$\left[x_1{,\mathrm{x}}_2\right]$	[2, 2]	[2, −2]	[0, 1]	/
$\left[\mathrm{u}\right]$	−0.5	1.5	−0.5	/

. In addition, the initial and the final conditions of the state and control variables are listed in the same table. In this case, only the initial condition of the state and control variables are given, with no limitation for the final condition.

To clearly illustrate the influence of uncertainty on the state trajectories, the optimized deterministic state trajectory and the state trajectories affected by uncertainty are shown in Figure 6 with different lines. On the basis of the uncertainty information of the design variables, MCS was employed to generate 100,000 samples, and the corresponding state trajectories were obtained using the trajectory deduction method, shown as green lines in the figure. The blue lines with circlular labels are the deterministic state trajectories.

As shown in Figure 6, the design variables and control strategy obtained by deterministic MDSDO did not satisfy the reliability requirements in Formula (14), which required the state variable $x_1$ to be larger than −0.4. Under the influence of uncertainty, the design variable could not remain constant, which makes the change in state trajectories and the constraint unsatisfactory. With the same illustration method as in Figure 6, the optimization result obtained by the proposed C-MDSDO approach is shown in Figure 7. The constraint violations were significantly improved, and only a small number of obtained state trajectories violated the condition, as shown in Figure 7. Moreover, the comparison of the state trajectories in Figure 8 shows that the proposed approach did not significantly change the control strategy but made only minor adjustments based on the original strategy.

To further demonstrate the advantages of the proposed method, the design variables scheme, the objective function performance, and the failure probability are compared in

Table 2. The optimization result comparison between MDSDO and C-MDSDO.

Method	Design Variable	Objective	Failure Probability
MDSDO	[0.0657, 1.0266]	1.9693	43.7%
RB-MDSDO [30]	[−0.008, 0.014]	1.9780	3.10%
C-MDSDO	[−0.0101, 0.0039]	1.9781	0.03%

. Herein, the optimization result of the RB-MDSDO approach, which was proposed in a previous study [30], is listed in the table. First, the proposed C-MDSDO method and the RB-MDSDO method both reduced the probability of failure. However, since the objective reliability index was 3 and the converted failure probability was approximately 0.27%, only the result of the C-MDSDO approach was deemed to be satisfactory, with a probability of 0.03%. Moreover, the failure probabilities of the original deterministic MDSDO approach and RB-MDSDO were 43.7% and 3.10%, respectively. Secondly, the performances of the objective function of the C-MDSDO and RB-MDSDO approaches were almost the same. Since the proposed approach has obvious advantages in terms of failure probability, it can be considered more reliable.

4.2. Glider Dynamic Soaring

This case is a modified optimization problem of the glider dynamic soaring. The optimization objective is the maximum $\overline{\rho }$, which is the fundamental parameter describing the performance of a wind gradient profile for a given glider [32]. The mathematical expression of the MDSDO problem for this case is shown in Formula (15), containing state variables $\left[\mathrm{x},y,h,V,\gamma ,\psi \right]$ and two control variables $\left[C_L,\mu \right]$.

$$\begin{array}{cc}\hfill \underset{\omega ,x,y,z,V,\gamma ,\psi ,C_L,\mu }{max}& \varphi =-\overline{\rho }=-\frac{\rho \mathrm{g}S}{2m{\omega }^2}\hfill \\ \hfill subjectto& \dot{x}=Vcos\left(\gamma \right)sin\left(\psi \right)+\omega Vsin\left(\gamma \right)\hfill \\ & \dot{y}=cos\left(\gamma \right)\hfill \\ & \dot{z}=Vsin\left(\gamma \right)\hfill \\ & \dot{\mathrm{V}}=-\frac{\rho S}{2m}\left(C_{D_0}+K{C}_L^2\right){\mathrm{V}}^2-gsin\left(\gamma \right)-\omega Vsin\left(\gamma \right)cos\left(\gamma \right)sin\left(\psi \right)\hfill \\ & \dot{\gamma }=\frac{\rho S}{2m}{C}_{L}Vcos\left(\mu \right)-\frac{1}{V}\left(gcos\left(\gamma \right)+\omega Vsin{\left(\gamma \right)}^{2}sin\left(\psi \right)\right)\hfill \\ & \dot{\psi }=\frac{1}{cos\left(\gamma \right)}\left(\frac{\rho S}{2m}{C}_{L}Vsin\left(\mu \right)-\omega sin\left(\gamma \right)cos\left(\psi \right)\right)\hfill \\ & {\mathrm{g}}_{\mathrm{i}}(·)\ge 0\mathrm{i}=1,2,3,4,5,6\hfill \\ & \frac{L}{mg}-5=\frac{\rho S}{2mg}{C}_L{V}^2-5\le 0\hfill \\ & X=\left[\mathrm{x},y,z,V,\gamma ,\psi \right],U=\left[C_L,\mu \right]\hfill \\ & {\omega }_{min}\le \omega \le {\omega }_{max}\hfill \\ & {\mathrm{X}}_{min}\le X\le X_{max}\hfill \\ & {\mathrm{U}}_{min}\le U\le U_{max}\hfill \end{array}$$

(15)

where $x$, $y$ are east and north positions, respectively; $z$ is the altitude, $V$ is the airspeed, $\gamma$ is the air-relative flight path angle, $\psi$ is the heading angle measured clockwise from the north, $C_L$ is the lift coefficient, and $\mu$ is the glider bank angle. In addition, ${\mathrm{g}}_{\mathrm{i}}$ represents the boundary conditions of state variables, $\omega$ is the mean of wind gradient slope, $C_{D_0}$ is the parasitic drag coefficient, $L$ is the lift force, $S$ is the windward area, and $m$ is the mass.

In the OCP, the value of the parameter $\omega$ is always constant and is calculated using Formula (16).

$$\omega =\frac{W_{\mathrm{x},max}}{h_{tr}}$$

(16)

where $W_{\mathrm{x},max}$ is the maximum horizontal wind, and $h_{tr}$ is the transition altitude at which the horizontal wind becomes constant.

In the practical situation, the parameter $\omega$ obeys the normal distribution [32]. Therefore, in reliability-based MDSDO research, this parameter is treated as the design variable and obeys the normal distribution $\omega \sim N\left({\mu }_{\omega },0.005\right)$. Moreover, the constraints of state variables are also transformed into reliability constraints, and the objective reliability index is ${\beta }_{obj}=2$.

The boundary conditions of all the variables are listed in

Table 3. The information of the design, state, and control variable of test 4.2.

Variable	Range		Initial	Final
	Max	Min
$\left[\omega \right]$	0.15	0.05	/	/
$\left[\mathrm{x},y,z,V,\gamma ,\psi \right]$	[1000, 1000, 750, 350, 1.31, 1.57]	[−1000, −1000, −20, 10, −1.31, −9.42]	[0, 0, 0, /, /, /]	[0, 0, 0, /, /, /]
$\left[C_L,\mu \right]$	[1.5, 1.31]	[−0.5, −1.31]	/	/

. Furthermore, the initial and the final condition of the state and control variables are listed in the same table. Angle-related parameters were calculated using the radian system (rad). In addition, the glider is required to return to the starting point within the specified time interval, and the time endpoint condition is $t_f\in \left[0,30\right]$.

Employing the same MCS-based presentation method as in the previous case, the state trajectories of deterministic MDSDO and the state trajectories of the design variables subject to uncertainty are shown in Figure 9. In the 1e5 MCS samples, 4967 samples did not satisfy all the constraints. Thus, the optimization result obtained using MDSDO did not meet the reliability requirements.

The state trajectories of the C-MDSDO optimization result are shown in Figure 10. As can be seen by comparing Figure 9 and Figure 10, the trajectory violation situation is significantly improved, especially in the trajectory of the state variable $x$ and the state variable $z$. The design variables scheme, the objective function performance, and the failure probability are compared in

Table 4. The optimization result comparison for MDSDO and C-MDSDO.

Method	Design Variable	Objective	Failure Probability
MDSDO	0.0637	−0.4461	49.67%
C-MDSDO	0.1293	−0.9050	0.09%

. Comparing the results of the MDSDO and the C-MDSDO approach, the failure probability dramatically decreased from 49.67% to 0.09%, which meets the reliability requirement. Furthermore, the optimization objective value decreased from −0.4461 to −0.9050.

As shown in Figure 11, the control strategy completely changed as a result of considering reliability, which made the glider turn back to the starting point in a shorter time, whereas the control curves were not as smooth as for the original solution. Additionally, as a result of the difference between the control strategy of the C-MDSDO and MDSDO approach, the trajectories of the state variables of the MDSDO method were different to those obtained by the C-MDSDO method.

4.3. Robot Arm

This case is derived from a point-to-point control problem of a three-degree-of-freedom robot arm [20], which contains six state variables $\left[{\theta }_1,{\theta }_2,{\theta }_3,{\omega }_1,{\omega }_2,{\omega }_3\right]$ and three control variables $\left[{\tau }_1,{\tau }_2,{\tau }_3\right]$. The structural diagram of this robot arm is shown in Figure 12.

In this paper, lengths of arms are treated as design variables $d=\left[L_1,L_2\right]$ for reliability-based MDSDO. The reliability constraint requires that the probability of the error between the actual final position of the robot arm and the target position is no more than 0.05 under the influence of uncertainty, and the converted reliability index is $\beta =3$. The mathematical expression of this optimization case is given in Formula (17).

$$\begin{array}{cc}\hfill \underset{d,u\left(t\right)}{min}& \varphi {=\mathrm{t}}_{\mathrm{f}}\hfill \\ \hfill subject to& \dot{x}=\mathrm{f}\left(d,x\left(t\right),u\left(t\right),t\right)\hfill \\ & x_0=\left[0,-1.5,0,0,0,0\right]\hfill \\ & x_{f_{MDSDO}}=\left[1,-1.95,1,0,0,0\right]\hfill \\ & prob\left[\Vert X-X_{obj}\Vert \ge 0.05\right]\le \Phi \left(-\beta \right)\hfill \\ & \beta =3,{\mathrm{L}}_1\sim N\left({\mu }_{L_1},0.006\right){,\mathrm{L}}_2\sim N\left({\mu }_{L_1},0.007\right)\hfill \\ & d=\left[L_1,L_2\right],x=\left[{\theta }_1,{\theta }_2,{\theta }_3,{\omega }_1,{\omega }_2,{\omega }_3\right],u=\left[{\tau }_1,{\tau }_2,{\tau }_3\right]\hfill \\ & d_{min}\le d\le d_{max}\hfill \\ & x_{min}\le x\le x_{max}\hfill \\ & u_{min}\le u\le u_{max}\hfill \end{array}$$

(17)

where $X$ represents the final position of the arm, and $X_{obj}$ represents the objective position, which is the final position of the deterministic MDSDO solution.

The boundary conditions of all the variables are listed in

Table 5. The information of the design, state, and control variable of test 4.3.

Variable	Range		Initial	Final
	Max	Min
$\left[L_1,L_2\right]$	[0.6, 1.2]	[0.3, 0.8]	/	/
$\left[{\theta }_1,{\theta }_2,{\theta }_3,{\omega }_1,{\omega }_2,{\omega }_3\right]$	[2.97, 5, 2.86, 5, 5, 10]	[−2.97, −5, −2.86, −5, −5, −6]	[0, −1.5, 0, 0, 0, 0]	[1, −1.95, 1, 0, 0, 0]
$\left[{\tau }_1,{\tau }_2,{\tau }_3\right]$	[7.5, 7.5, 7.5]	[−7.5, −7.5, −7.5]	/	/

. In addition, the initial and the final conditions of the state and control variables are listed in the same table. Angle-related parameters were calculated using the radian system (rad). In addition, the arm was required to reach the target point within the specified time interval, and the time endpoint condition was $t_f\in \left[0,0.5\right]$.

In order to obtain the target position of the robot arm, it was necessary to convert the relative angular coordinates in Table 5 into the coordinates in the Cartesian coordinate system. The original arm length [0.5, 0.98] was substituted into the coordinate conversion (18) to obtain the target end position $X_{t_f}$. Furthermore, during reliability-based MDSDO, the end position of the robot arm also needed to be converted and compared with $X_{t_f}$, rather than directly using the relative angle $x_{f_{MDSDO}}=\left[1,-1.95,1\right]$.

$$\{\begin{array}{l}{X}_{t_f}{}_x{=L}_{1}cos\left({\theta }_1\right)+L_{2}cos\left({\theta }_1+{\theta }_2\right)\\ X_{t_f}{}_y{=L}_{1}sin\left({\theta }_1\right)+L_{2}sin\left({\theta }_1+{\theta }_2\right)\end{array}$$

(18)

The same state trajectory presentation method was used here to demonstrate the effect of uncertainty on the state trajectory of the robot arm in Figure 13. With the influence of uncertainty, the angle of the final state of the robot arm varied greatly. Moreover, the state variables $\left[{\omega }_1,{\omega }_2,{\omega }_3\right]$ represent the angular acceleration of the arm, which means that the arm was not stationary at the endpoint under nondeterministic conditions. When the end position was converted using Formula (18), approximately 30.8% of the MCS samples indicated that the end position error was greater than the threshold, which demonstrates the unreliability of the MDSDO solution and shows that reliability-based MDSDO is required.

As in the previous case, the trajectories of the state variables of the solution obtained using the C-MDSDO approach are shown in Figure 14. The Figure 14 shows that the differences between the deduced state trajectories of each MCS sample are reduced. Moreover, the control strategy of the C-MDSDO solution became easier to achieve in practice, as shown in Figure 15. In contrast, the strategy from the original MDSDO solution exhibited high levels of fluctuation, which means that the control performance of the controller was exceptionally demanding.

The design variables scheme, the objective function performance, and the failure probability are compared in

Table 6. The optimization result comparison for MDSDO and C-MDSDO.

Method	Design Variable	Objective	Failure Probability
MDSDO	[0.3000, 0.8000]	0.4073	30.08%
C-MDSDO	[0.4483, 0.7694]	0.3422	0.13%

. Using Formula (18) to transform the coordinates of the end position of the MCS random samples generated according to the solution obtained by C-MDSDO, only 0.13% of the sample end position errors were greater than the threshold, while the probability of the original MDSDO solution was 30.08%. Moreover, the performance of the objective function decreased from 0.4073 (MDSDO) to 0.3422 (C-MDSDO). It is worth noting that the solution using the original MDSDO approach was inferior to that of C-MDSDO, not only in terms of the performance of the objective function but also in terms of the control strategy. The reason for this is that the reliability constraint in this example is a range constraint, which causes the design region of C-MDSDO to increase in size as compared to that of the MDSDO condition. Another target point of the robot arm was found in the feasible domain by C-MDSDO, which means the reliability meets the objective reliability constraint.

4.4. Discussion

All the test results met the set reliability index; however, the performance of the objective function was different. In Case 4.1, the performance value decreased slightly from 1.9693 to 1.9781. In the other two cases, the variation was different, and this is related to the characteristics of the optimization problem. In Case 4.2 and Case 4.3, the state trajectories and control strategies of C-MDSDO were significantly different from those of MDSDO, which causes a significant variation in the performance of the objective function. Comparing the termination conditions of the three cases, we found that neither Case 4.2 nor Case 4.3 had a specific termination time point as compared to Case 4.1. Therefore, in the last two cases, the algorithm was able to obtain a solution that satisfied the reliability constraint only by applying a more aggressive control strategy. However, as mentioned previously, the range reliability constraint of Case 4.3 increased the size of the feasible domain of C-MDSDO as compared to that of MDSDO, which improved the objective performance but shortened the control time. It is worth mentioning that, in Case 4.1, the RB-MDSDO method did not strictly meet the required reliability index, but C-MDSDO succeeded, as compared with the limited available research. Moreover, the performance of the objective function for C-MDSDO and RB-MDSDO was almost equal.

In order to test the robustness of the proposed approach, several repeated experiments were performed for each case. In

Table 7. Repeated experiment results of cases (number of experiments = 30).

Method	Case 4.1	Case 4.2	Case 4.3
Objective performance	1.9781 ± 1.0259 × 10−8	−0.9050 ± 1.3696 × 10−32	0.3422 ± 3.9248 × 10−6
Failure probability	0.0329% ± 5.3498 × 10−5%	0.0937% ± 7.6118 × 10−10%	0.1339% ± 1.6242 × 10−4%

, the objective function and failure probability of each case with multiple experiments are presented in the form of ”mean ± variance”. In the table, the variance of the objective performance and failure probability is small for each case, which indicates that the proposed approach converges stably to the same solution in repeated experiments. The data in the table demonstrate that the proposed approach can consistently obtain the design that meets the reliability requirement.

5. Conclusions

In this paper, the C-MDSDO method is proposed to solve the reliability-based MDSDO problem. Several cases were employed to test the effectiveness of the proposed method, and the results from relevant studies were also used for comparison. From the test results, the proposed method was shown to be efficient and reliable. Specifically, in the numerical case, i.e., Case 4.1, when compared with the RB-MDSDO approach, the optimization results of C-MDSDO were more reliable, while the performances of the objective function were similar. In the engineering cases, i.e., Case 4.2 and Case 4.3, faster control strategies were obtained to meet the reliability constraint functions, but the performance of the objective function and the smoothness of the control strategy degraded. In addition, the state trajectories of the optimization objects changed significantly as compared with the MDSDO solution.

In summary, the main contributions of this paper are as follows: (1) an improved sequential approach, C-MDSDO, is proposed to solve the problem of the unexpected failure of the MDSDO solution due to design variable uncertainty; (2) an optimization problem was designed based on the crossrate of the objective reliability index sample to rapidly obtain the shifting vector, which is employed to shift the solution into the reliability domain.

In addition, the approach proposed in this paper is only applicable to cases in which the differential equation of the dynamic discipline is known. For cases in which the differential equation information is unknown, the surrogate model technique is required to form the black-box-based approach. In future research, the influence of uncertainty on the control strategy needs to be further investigated to enhance the generalizability of the reliability-based MDSDO theory.