A Novel Surrogate Model-Based Solving Framework for the Black-Box Dynamic Co-Design and Optimization Problem in the Dynamic System

Abstract

When encountering the black-box dynamic co-design and optimization (BDCDO) problem in the multidisciplinary dynamic system, the finite difference technique is inefficient or even infeasible to provide approximate numerical gradient information for the optimization algorithm since it requires numerous original expensive evaluations. Therefore, a solving framework based on the surrogate model of the state equation is introduced to optimize BDCDO. To efficiently construct the surrogate model, a sequential sampling method is presented on the basis of the successive relative improvement ratio. Meanwhile, a termination criterion is suggested to quantify the convergence of the solution. Ultimately, the newly proposed sampling strategy and termination criterion are incorporated into the BDCDO solving framework to optimize two numerical examples and two engineering examples. The results demonstrate that the framework integrating the proposed sampling strategy and termination criterion has the best performance in terms of the accuracy, efficiency, and computational budget compared to the existing methods.

1. Introduction

Differing from the dynamic optimization problem (DOP), also known as the optimal control problem, in which only the control strategy decision is optimized to improve the performance of the dynamic system [1,2,3], the dynamic co-design and optimization (DCDO) problem accounts for the bi-directional dependency of physical system design and control system design and includes two types of design variables: plant (or physical) and control [4]. Two categories of co-design methods, nested (or multi-layer optimization) method and simultaneous method, are developed and deployed on the DCDO problem in the engineering applications [5,6,7,8]. In those two co-design methods, though the optimization structures are different, the DCDO problem should be transcribed into a finite-dimensional nonlinear programming (NLP) problem via direct transcription [9] to optimize.

However, when solving the DCDO problem of the sophisticated dynamic systems involving multiple disciplines or multiple subsystems, co-design schemes inevitably encounter obstructions such as high computational consumption incurred from the time-consuming system simulations [10]. Moreover, in some engineering practices, dynamic system models are constructed by industrial simulation software or platforms, and the explicit equations of the state in dynamic systems expressed by differential algebraic equations cannot be extracted directly from the dynamic models [11]. The DCDO involving such dynamic system models is referred to as the black-box dynamic co-design and optimization (BDCDO) problem. Due to the lack of explicit state equation, the finite-difference technique, rather than the automatic differentiation technique, provides approximate Jacobian information matrices for NLP solvers to iteratively optimize BDCDO, which demands significant computational resources to evaluate the dynamic systems and definitely further increases the computational budget of the co-design schemes.

The response surface methods (RSMs) have been proven to be effective tools to address computationally expensive problems in complex static black-box systems [12,13,14,15,16]. Therefore, RSMs are introduced to approximate the black-box dynamic system to alleviate the number of the original system valuations and save computational costs. In attempting to reduce the modeling difficulty and preserve the dynamic properties of the dynamic system, Deshmukh et al. [17] presented the derivative function surrogate modeling methodology to construct surrogate models for the derivative functions of the dynamic system rather than construct surrogate models of the whole system responses [18,19,20]. In order to build high-fidelity model surrogate models of the derivative functions, Deshmukh et al. [17] used the Latin Hypercube Sampling (LHS) method for sequential sampling in the minimal hypercube space containing the current optimal trajectory to update the surrogate model. Lefebvre et al. [21] averaged the errors of the KRG model after integrating the errors along the current optimal trajectory and then used the values of state variables obtained by inverse error integral of each segment as new samples to update the KRG model. Qiao et al. [22] proposed EFDC sampling method based on KRG model to filter the current trajectory discrete points after error analysis, combine the spatial distance to cluster these discrete points, and select the points with the largest prediction error to update the model. At the same time, some reasonable solution termination criteria also have been investigated to avoid redundant iterations, which contribute little to the accuracy improvement of the solution result. Deshmukh et al. [17] determined whether the solution process stops or not based on the discrepancy between the current and previous iterates. Lefebvre et al. [21] proposed a new metric, dynamical mismatch, to identify whether the solution process is terminated or not. However, computing the dynamical mismatch needs additional sample points and consumes more computational resources. Qiao et al. [22] adopted the accuracy of the surrogate model and the successive relative improvement of the objective function as the stopping plan to assess the convergence of the solution process more comprehensively.

Admittedly, the BDCDO solving framework combined with the above-mentioned sampling strategies and termination criteria indeed reduce the number of samples for constructing the surrogate models of derivative functions to different degrees, but the efficiency of modeling and the robustness of the optimal solution still require further improvement. To this end, a new sequential sampling strategy based on the successive relative improvement ratio of the discrete trajectory points and maximizing distances between the nearest sample points, called SRIRMD, is proposed in this work to effectively improve the accuracy of the surrogate models by selecting the points with a large successive relative improvement ratio among the discrete trajectory points as new samples. Meanwhile, maximizing the minimum distances between new samples and existing samples can ensure the uniform distribution of all sample points. In addition, according to the fundamental observation that the state trajectories tend to coincide during the solving process of dynamic optimization problems, a new termination criterion, named the state trajectory overlap ratio (STOR), is presented to quantify the convergence and intuitively reflect the convergence trend of the solution in the iterative process. Finally, two numerical examples, one 3-DOF robot co-design and optimization problem and one horizontal axis wind turbine co-design and optimization problem, are solved by the means of the BDCDO solving framework integrated with the SRIRMD sampling strategy and the STOR termination criterion. The results demonstrate that the BDCDO solving framework combining the SRIRMD sampling strategy and the STOR termination criterion has the best performance compared to existing methods and can obtain more accurate and robust solutions with fewer sample points, improving the solution efficiency and reducing the computational budget.

The rest of this paper is organized as follows. Section 2 reviews the dynamic optimization problem and its direct solving method and the Kriging technique. Section 3 introduces the BDCDO solving framework combining the new sampling strategy and termination criterion. Section 4 verifies the feasibility and efficiency of the BDCDO solving framework through two numerical examples and two engineering examples. The conclusions are revealed in the last section.

2. Background

2.1. Dynamic Co-Design and Optimization Problem and Its Direct Solving Method

The objective of the dynamic co-design and optimization (DCDO) problem is to solve the optimal input vectors of physical design parameters $x_p^{*}$ and control variables $u^{*}(t)$ that minimizes the system performance index. The standard form of DCDO based on the simultaneous method is described as follows [23,24]:

$$\begin{array}{ll}\underset{x_p,u(t)}{\min }& J(x_p,u(t))=\varphi (x_p,\xi (t_0),\xi (t_f),t_0,t_f)+{\displaystyle {\int }_{t_0}^{t_f}L(x_p,\xi (t),u(t),t)dt}\\ s.t.& \dot{\xi }(t)=f(x_p,\xi (t),u(t),t)\\ & \Psi (\xi (t_0),t_0,\xi (t_f),t_f)=0\\ & g(x_p,\xi (t),u(t),t)\le 0\\ & u_L\le u(t)\le u_U\\ & x_L\le x_p\le x_U\end{array}$$

(1)

where $J$ is the response of a cost function that consists of a Mayer term $\varphi (\cdot )$ and Lagrange term $L(\cdot )$. $t\in [t_0,t_f]$ denotes the time horizon, and $t_0$ and $t_f$ are initial time and terminal times. $\xi (t_0)$ and $\xi (t_f)$ indicate the initial and terminal states of the system, while $\xi (t)$ and $u(t)$ are vectors of the state variables and control inputs at moment $t$. DCDO is subject to several different constraints such as state equations $\dot{\xi }(t)=f(x_p,\xi (t),u(t),t)$, boundary constraints $\Psi$, and path constraints $g$. Among these three constraints, the state equations and path constraints are continuous constraints that must be satisfied during the entire time period, while the boundary constraints are discrete constraints that only need to be satisfied at the moments $t_0$ and $t_f$. Finally, the vectors of the physical design parameters $x_p$ and control inputs $u(t)$ enforce the interval constraints $[x_L,x_U]$ and $[u_L,u_U]$, respectively.

In the direct solving approach of the DCDO, the vectors of the state variables $\xi (t)$, plant design parameters $x_p$, and control inputs $u(t)$ are discretized at the time grid nodes by means of the DT. Thus, the DCDO is transformed into an NLP, as shown below.

$$\begin{array}{ll}\underset{x_p,\Xi ,\Theta }{\min }& J(x_p,\Xi ,\Theta )=\varphi (x_p,\xi (t_0),\xi (t_f),t_0,t_f)+{\displaystyle \sum _{k=0}^N{\omega }_k}\cdot L(t,x_p,\Xi (k),\Theta (k))\\ s.t.& D\cdot \Xi =f(x_p,\Xi ,\Theta )\\ & \Psi (\xi (t_0),t_0,\xi (t_f),t_f)=0\\ & g(x_p,\Xi ,\Theta )\le 0\end{array}$$

(2)

where $\Xi$ and $\Theta$ are the discrete matrices of the state variables and control inputs in the time domain, ${\omega }_k$ is the integration weight, and $D$ is the differential matrix specified in different pseudospectral methods [25,26]. When solving the BDCDO of the dynamic system based on the surrogate models of the derivative functions, the right hand-side functions of the state equations are approximated by the surrogate models. Hence, replacing the derivative function $f(\cdot )$ in Equation (2) with the surrogate model $\widehat{f}(\cdot )$ forms a computationally inexpensive expression:

$$D\cdot \Xi =\widehat{f}(x_p,\Xi ,\Theta )$$

(3)

2.2. Kriging Technique

In view of the efficiency and effectiveness of the Kriging technique in approximating low-dimensional problems, as well as providing prediction errors at arbitrary points [27], the Kriging technique is adopted to construct the surrogate model in the BDCDO. When training the model, the state variables, physical design parameters, and control variables $X={[x_p,\xi ,u]}^T$ are fed into the model, and the output are the valuations of the derivative functions $f(X)$ in the state equation.

The standard Kriging model has two main components: the regression function and the stochastic process [28]. The expression is formulated as follows:

$$\widehat{f}(X)=F^T\left(X\right)\mathsf{\beta }+z(X)$$

(4)

where $F^T(X)$ contains a series of regression functions, and $\mathsf{\beta }$ is a trend coefficient vector. $z(X)$ denotes a random process with a mean of 0 and variance ${\sigma }_z^2$. The spatial covariance function between the stochastic processes $z(X_i)$ and $z(X_j)$ can be expressed as

$$\mathrm{cov}[z(X_i),z(X_j)]={\sigma }_z^{2}G(\mathsf{\theta },X_i,X_j)$$

(5)

where $X_i$ and $X_j$ are two different points in the design space, $G$ is the Gaussian spatial correlation function, and $\mathsf{\theta }$ is the corresponding vector of correlation coefficients. The values of $G$ and $F$ at sample points constitute matrices $R=G(X_i,X_j)$ and $S=F(X_i)$. According to the unbiased estimator theory, the least squares solution of the regression model is

$$\widehat{\beta }={(S^T{R}^{-1}S)}^{-1}{S}^T{R}^{-1}Y$$

(6)

and the maximum likelihood estimation of variance is

$${\widehat{\sigma }}_z^2=\frac{1}{m}{(Y-S\widehat{\beta })}^T{R}^{-1}(Y-S\widehat{\beta })$$

(7)

where $Y=[\widehat{f}(X_i)], i=1,2,\dots ,m$.

Thus, Equation (4) is translated into

$$\widehat{f}(X)=F^T\left(X\right)\beta +r^T\left(X\right)R^{-1}(Y-S\widehat{\beta }), r\left(X\right)={[G(\theta ,X,X_1),\dots ,G(\theta ,X,X_m)]}^T$$

(8)

Suppose $d=S^T{R}^{-1}r-S$, the predictor of the mean square error (MSE) is calculated as follows:

$$MSE={\sigma }_z^2(1+d^T{(S^T{R}^{-1}S)}^{-1}d-r{R}^{-1}r)$$

(9)

3. The BDCDO Solving Framework

The BDCDO solving framework based on the surrogate model of the derivative function mainly consists of two parts: refining the surrogate model for the right hand-side function of the state equation and solving the black-box dynamic co-design and optimization (BDCDO) problem based on the surrogate model. In the loop of approximating the derivative function, the initial model is built according to the initial samples set and then updated by the sequential sampling method. In the loop of solving the BDCDO, the BDCDO is discretized into NLP at time grid nodes, then the approximate Jacobian information matrices based on the model are delivered for the SQP algorithm to solve the NLP. For the purpose of efficiently constructing the surrogate model of the derivative function, the new sampling method based on the successive relative improvement ratio of the discrete trajectory points and maximize the distances between the sample points, termed SRIRMD, is elaborated in this section. Meanwhile, to quantify the convergence and intuitively reflect the convergence trend of the solution, a new termination criterion on the basis of the fact that the state trajectories tends to coincide during the solving process, called the state trajectory overlap ratio (STOR), is also introduced in this section. Finally, the newly proposed sampling strategy, SRIRMD, and termination criterion, STOR, are integrated into the solving framework of BDCDO.

3.1. Adaptive Sequential Sampling Strategy

In the static optimization problem, the optimal result is a point; hence, single sampling strategies that focus on improving the accuracy of the surrogate model near the current optimal point are widely used, while batch sampling strategies are evidenced to be unable to enhance the performance of updating model [29]. However, the scenario is quite different in the dynamic optimization problem since the optimal results are time-dependent control curves and the corresponding state trajectories. Adding a new sample point during each iteration contributes less to improving the convergence rate of the state trajectory but instead increases the number of modeling and optimization throughout the solution process and consumes more computational resources.

Motivated by the sequential sampling strategy in the static optimization problem, the adaptive sampling strategy applied to the DOP should concentrate on how to improve the accuracy of the local region where the current state trajectory is located in order to avoid redundant sampling in the uninteresting area. To this end, it is necessary to pick informative points in the vicinity of the current trajectory to update the surrogate model of the derivative function. Fortunately, the discrete trajectory points (DTPs) on the current state trajectory offer a large number of candidate samples. However, it is not wise to add all DTPs to the samples set, because (a) the expensive evaluation of all DTPs requires a lot of computational effort, and (b) closer DTPs are less helpful to improving the accuracy but increase the complexity of the model. Hence, a new sequential sampling strategy based on the successive relative improvement ratios of the DTPs and maximizing the distances between the sample points, called SRIRMD, is presented in this work. The SRIRMD sampling strategy prioritizes picking the points with large successive relative improvement ratios among the DTPs as new samples to refine the surrogate model. Meanwhile, maximizing the minimum distances between new samples and existing samples can guarantee that all sample points are distributed uniformly. The specific steps of the SRIRMD method are listed below.

Step 1: Obtain the initial values $V_I$ and the optimal values $V_O$ at the discrete points of the current state trajectory and calculate the successive relative improvement ratios (srir) of all discrete points.

$$SRIR=\left\{srir|srir=\left|V_O^i-V_I^i\right|/\left|V_O^i\right|,i=1,2,\dots ,m\right\}$$

(10)

where $m$ is the number of discrete points of the current trajectory.

Step 2: Generate the candidate samples set $S_{SRIR}$ by eliminating the points whose srir is smaller than the allowable deviation factor $\beta$ from the DTPs.

$$S_{SRIR}=\left\{DTPi\right|sriri\ge \beta ,sriri\in SRIR\}$$

(11)

Step 3: Calculate the distance matrix $D_{\min }$ of all candidate samples in $S_{SRIR}$ to the nearest point in the existing samples set $S^{(l)}$.

Step 4: Select new sample $x_{new}$ by the following sampling criterion.

$$\max srir(x_{new})\cdot D\min (x_{new})$$

(12)

Step 5: Add the new sample $x_{new}$ into samples set $S^{(l)}$ and delete $x_{new}$ in $S_{SRIR}$.

Step 6: Repeat Step 3–5 until the required number of samples are picked, and the samples set is updated to $S^{(l+1)}$.

Step 7: Terminate the current round of sampling.

Figure 1 exhibits the process of selecting new sample points from the DTPs employing the SRIRMD strategy. The dotted lines are the previous state trajectories, the dots are the discrete points of those previous state trajectories, and the red solid lines are the current optimal state trajectories. The new points sampled by SRIRMD method are plotted in the form of black stars, and these samples are mainly distributed in the areas with large successive relative improvement ratios while being evenly spread.

3.2. Termination Criterion: The State Trajectory Overlap Ratio

The state trajectory overlap ratio, STOR, is designed to intuitively assess the convergence of state trajectories and serve as the termination criterion. The formula for STOR is expressed as follows:

$$\mathrm{A}={\displaystyle \prod _{i=1}^d{\alpha }_i}$$

(13)

where $d$ is dimension of state variables. ${\alpha }_i$ means the state component trajectory overlap ratio, which reflects the trajectory overlap ratio of the $i_{th}$ state component in two successive iterations. As shown in Figure 2, $S^{*}$ is the trajectory of the state component ${\xi }_i$ obtained in the previous iteration, and $S^{**}$ is the trajectory of ${\xi }_i$ in the current iteration. The ${\alpha }_i$ of the state component ${\xi }_i$ can be calculated by the lengths of $\stackrel{⌢}{\mathrm{A}B}$ and $\stackrel{⌢}{\mathrm{A}C}$, and the specific formula is expressed in Equation (14).

$${\alpha }_i=\frac{L_{\stackrel{⌢}{\mathrm{A}B}}}{L_{\stackrel{⌢}{\mathrm{A}C}}}$$

(14)

Equation (14) denotes the formula for calculating the ${\alpha }_i$ in the state space. However, the lengths of the trajectories $\stackrel{⌢}{\mathrm{A}B}$ and $\stackrel{⌢}{\mathrm{A}C}$ are not convenient to calculate in the actual computation procedure. To facilitate the calculation, Equation (14) can be mapped from the state space to the time domain since the state trajectories and time points correspond (i.e., each time point corresponds to a specific trajectory value). Hence, Equation (14) is replaced by Equation (15) to compute ${\alpha }_i$.

$${\alpha }_i=\frac{L_{\stackrel{⌢}{\mathrm{A}B}}}{L_{\stackrel{⌢}{\mathrm{A}C}}}$$

(15)

where $T_0$ and $T_1$ are represent the overlap and non-overlap periods of trajectories $S^{*}$ and $S^{**}$ in the time domain $T$, respectively. $T_0$ and $T_1$ can be computed by the following expression.

$$\begin{array}{l}{T}_0=\left\{t|\frac{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t),u^{**}(t);S{M}^{**})-{\dot{\xi }}_i(x_p^{*}(t),{\xi }^{*}(t),u^{*}(t);S{M}^{*})\right|}{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t),u^{**}(t);S{M}^{**})\right|}>\beta ,t\in T\right\}\\ T_1=\left\{t|\frac{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t),u^{**}(t);S{M}^{**})-{\dot{\xi }}_i(x_p^{*}(t),{\xi }^{*}(t),u^{*}(t);S{M}^{*})\right|}{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t),u^{**}(t);S{M}^{**})\right|}\le \beta ,t\in T\right\}\end{array}$$

(16)

where $S{M}^{*}$ and $S{M}^{**}$ are surrogate models of derivative function in the previous and current iterations. $\beta$ is the allowable deviation factor of the trajectory, which indicates trajectories $S^{*}$ and $S^{**}$ are also regarded as coincide when their error rate at moment $t$ are not greater than $\beta$. It is worth noting that $\beta =0.01$ in this research. The time domain $T$ is a continuous time series, and it could be uniformly discretized as $T=\{t_0,t_1,\dots ,t_{\tau },\dots ,t_f\}$ to accelerate computations. Thus, Equation (16) is converted into Equation (17).

$$\begin{array}{l}{T^{\prime }}_0=\left\{t_{\tau }|\frac{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t_{\tau }),u^{**}(t_{\tau });S{M}^{**})-{\dot{\xi }}_i(x_p^{*}(t),{\xi }^{*}(t_{\tau }),u^{*}(t_{\tau });S{M}^{*})\right|}{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t_{\tau }),u^{**}(t_{\tau });S{M}^{**})\right|}>\beta ,\tau =0,1,2,\dots ,f\right\}\\ {T^{\prime }}_1=\left\{t_{\tau }|\frac{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t_{\tau }),u^{**}(t_{\tau });S{M}^{**})-{\dot{\xi }}_i(x_p^{*}(t),{\xi }^{*}(t_{\tau }),u^{*}(t_{\tau });S{M}^{*})\right|}{\left|{\dot{\xi }}_i(x_p^{**}(t),{\xi }^{**}(t_{\tau }),u^{**}(t_{\tau });S{M}^{**})\right|}\le \beta ,\tau =0,1,2,\dots ,f\right\}\end{array}$$

(17)

Furthermore, Equation (14) for calculating ${\alpha }_i$ is transformed into Equation (18).

$${\alpha }_i=\frac{n_{T_1^{\prime }}}{n_{T_0^{\prime }}+n_{T_1^{\prime }}}$$

(18)

where $n_{T_0^{\prime }}$ and $n_{T_1^{\prime }}$ denote the number of elements in the sets $T_0^{\prime }$ and $T_1^{\prime }$, respectively.

According to the above series of transformations, including mapping and discretization, the calculation of ${\alpha }_i$ is finally transformed into the statistics of the elements in the sets. Obviously, the denser the time domain $T$ is divided in Equation (17), the more accurate ${\alpha }_i$ is in Equation (18). At the same time, the following two theorems are derived on the basis of the above definition about $\mathrm{A}$ and ${\alpha }_i$.

Theorem 1.

The sufficient condition for the convergence of the state trajectory overlap ratio$\mathrm{A}$is that all the state component trajectory overlap ratio${\alpha }_i$converge.

Proof: 

From the definitions of $\mathrm{A}$ and ${\alpha }_i$, it follows that $\mathrm{A}\in [0,1],{\alpha }_i\in [0,1]$. From the perspective of the mathematics, convergences of $\mathrm{A}$ and ${\alpha }_i$ implies that the values of $\mathrm{A}$ and ${\alpha }_i$ converge to 1, i.e.,

$$\mathrm{A} and {\alpha }_i converge<=>\mathrm{A}\to 1,{\alpha }_i\to 1$$

(19)

• If the dimension of the state variables $d=1$, $\mathrm{A}={\alpha }_1$, the theorem is established.
• If the dimension of the state variables $d>1$, suppose

$${\mathrm{A}}^k={\displaystyle \prod _{i=1,2,..,k-1,k+1,..,d}{\alpha }_i}$$

(20)

When $k=d$,

$$\mathrm{A}={\mathrm{A}}^{d-1}\cdot {\alpha }_d$$

(21)

where ${\mathrm{A}}^{d-1}\in [0,1]$ and ${\alpha }_d\in [0,1]$. Without considering the specific value of ${\mathrm{A}}^{d-1}$, only ${\alpha }_d\to 1$, then $\mathrm{A}$ may converge to 1. Conversely, if ${\alpha }_d\bcancel{\to }1$, then $\mathrm{A}$ definitely does not converge to 1 whether or not ${\mathrm{A}}^{d-1}$ converges to 1. Similarly, let k be $d-2,d-3,..,2$, respectively, so that the sufficient condition for $\mathrm{A}\to 1$ is ${\alpha }_i\to 1,i=1,2,..,d$. That means the sufficient condition for the convergence of $\mathrm{A}$ is that all ${\alpha }_i$ converge. □

Theorem 2.

When the state trajectory overlap ratio$\mathrm{A}$converges, every state component trajectory overlap ratio${\alpha }_i$is greater than or equal to $\mathrm{A}$.

Proof: 

• If the dimension of the state variables $d=1$, $\mathrm{A}={\alpha }_1$, the theorem is established.
• If the dimension of the state variables $d>1$, Equation (18) can be obtained based on Equations (20) and (21).

$${\alpha }_k=\frac{\mathrm{A}}{{\mathrm{A}}^k}$$

(22)

Suppose $\mathrm{A}$ is converged, then ${\mathrm{A}}^k$ is also converged (i.e., ${\mathrm{A}}^k\to 1$), since ${\mathrm{A}}^k\subset \mathrm{A}$. Obviously, ${\alpha }_k\ge \mathrm{A}$ can be deduced from Equation (22). Similarly, let k be $1,2,\dots ,d$, respectively, then ${\alpha }_i\ge \mathrm{A},i=1,2,\dots ,d$ can be proved. That means when $\mathrm{A}$ converges, all ${\alpha }_i$ are greater than or equal to $\mathrm{A}$.

According to the definition of the state trajectory overlap ratio $\mathrm{A}$ and the proving procedures of the above two theorems, it could be concluded that taking $\mathrm{A}$ as the convergence criterion for the BDCDO solving framework has the following advantages: (a) directly reflects the convergence trend of state trajectory in the iterative solving process, and $\mathrm{A}\in \left[0,1\right]$, which also quantifies the convergence; and (b) guarantees the convergences of state component trajectories. The state solution of the BDCDO converges only if all the state component trajectories converge. □

3.3. The BDCDO Solving Framework Combined with SRIRMD and STOR

By integrating the SRIRMD sampling strategy and the STOR termination criterion proposed in this work into the BDCDO solving framework, a new BDCDO solving method, named SRIRMD-STOR, is generated as depicted in Figure 3. The input and output of the SRIRMD-STOR method are shown in

Table 1. Input and output of the SRIRMD-STOR method.

Input

The upper and lower bounds of physical design parameters $x_p$, state variables $\xi$ and control inputs $u$. The initial guess values $[x_p^{(0)},{\Xi }^{(0)},{\Theta }^{(0)}]$ of physical design parameters $x_p$, state variables $\xi$ and control inputs $u$. The termination criterion threshold ${\mathrm{A}}_0$ for the SRIRMD-STOR method. The solving tolerance and max iteration for the DOP solvers [30].

Output

The optimal design point of physical design parameters $x_p$. The optimal trajectories of state variables $\xi$. The optimal curves of control inputs $u$.

, and the specific realization of the SRIRMD-STOR method is listed in

Table 2. The specific realization of the SRIRMD-STOR method.

The SRIRMD-STOR Method

Step 1: Apply LHS method to sample initial points set $S^0$ in the design domain consisting of feasible regions of physical design parameters $x_p$, state variables $\xi$ and control inputs $u$. Step 2: Construct the initial surrogate model ${\widehat{f}}^{(0)}$ of the derivative function by Kriging technique [31] with the initial samples set $S^0$. Step 3: Transcribe the BDCDO into the NLP at the time grid nodes via DT, then solve NLP based on the initial guess values of $[x_p^{(0)},{\Xi }^{(0)},{\Theta }^{(0)},{\widehat{f}}^{(0)}]$ and obtain the current optimal plant design parameters, state trajectories, control curves, and performance index $[x_p^{(1)},{\Xi }^{(1)},{\Theta }^{(1)},J^{(1)}]$. Step 4: Calculate the state component trajectory overlap ratios ${\alpha }_i$ of all state variables according to the initial guess trajectories and the current optimal trajectories, then calculate the state trajectory overlap ratio $\mathrm{A}$. If $\mathrm{A}>{\mathrm{A}}_0$, terminate the solving process; otherwise, go to Step 5. Step 5: Employ the SRIRMD strategy to select new samples $x_{new}$ from the current DTPs, update the samples set $S^1$, and rebuild the surrogate model ${\widehat{f}}^{(1)}$. Step 6: Update the time grid nodes using the grid optimization algorithm and translate the BDCDO into the NLP at the new time grid nodes. Step 7: Solve NLP based on the current values of plant design parameters; state trajectories and control inputs, and current model $[x_p^{(l)},{\Xi }^{(l)},{\Theta }^{(l)},{\widehat{f}}^{(l)}]$; and acquire the latest optimal plant parameters, state trajectories, control curves, and performance index. Note: $l$ starts from 1. Step 8: Calculate ${\alpha }_i$ and $\mathrm{A}$. If $\mathrm{A}>{\mathrm{A}}_0$, stop the solving process; otherwise, go to Step 5.

.

4. Numerical Examples

In this section, two numerical examples are used to verify the feasibility and effectiveness of the BDCDO solving framework combined with the SRIRMD sampling strategy and the STOR termination criterion, also called the SRIRMD-STOR method for short. The first numerical example is a classical multidimensional nonlinear dynamic optimization problem involving six state variables and two control inputs. The second example is a numerical nonlinear dynamic co-design and optimization problem that is sensitive to the accuracy of the grid optimization, i.e., it has high requirements on the number and location of the discrete grid nodes of the optimization problem. To solve those two numerical tests, STOR $\mathrm{A}$ is regarded as the termination criterion, and the threshold of $\mathrm{A}$ is set to 0.95.

4.1. Example 1: A Mathematical Nonlinear Dynamic Optimization Problem

The mathematical model of the this nonlinear dynamic optimization problem [32] can be described as follows:

$$\begin{array}{ll}\min & J={\displaystyle {\int }_0^{t_f}[L_2+I_2+\frac{1}{2}{B}_1{u}_1^2+\frac{1}{2}{B}_2{u}_2^2]dt}\\ s.t.& \dot{S}=\mathrm{\Lambda }-\frac{13}{30000}S{I}_1-\frac{0.029}{30000}S{I}_2-0.0143S\\ & \dot{T}=2u_1{L}_1-0.0143T+(1-0.5(1-u_2))I_1-\frac{13}{30000}T{I}_1-\frac{0.029}{30000}T{I}_2\\ & {\dot{L}}_1=\frac{13}{30000}S{I}_1-0.6143L_1-2u_1{L}_1+0.4(1-u_2)I_1+\frac{13}{30000}T{I}_1-\frac{0.029}{30000}{L}_1{I}_2\\ & {\dot{L}}_2=0.1(1-u_2)I_1-1.0143L_2+\frac{0.029}{30000}(S+L_1+T)I_2\\ & {\dot{I}}_1=0.5L_1-1.0143L_1\\ & {\dot{I}}_2=L_2-0.0143I_2\\ & N=S+L_1+L_2+I_1+I_2+T=30000\end{array}$$

(23)

where the state variables $\xi ={[S,T,L_1,L_2,I_1,I_2]}^T$. The upper and lower bounds of control inputs $u={[u_1,u_2]}^T$ are $u_U={[0.95,0.95]}^T$ and $u_L={[0.05,0.05]}^T$. The finial time $t_f$ is set as 5. To solve this DOP based on the mathematical model, the maximum number of iterations of the NLP solver is set to 15 and the solution accuracy is set to ${10}^{-4}$ to obtain the exact solution 5152.1. Although the mathematical model provides the explicit expressions of the state equation, it can still be treated as a black-box dynamic optimization problem and solved using the BDCDO solving framework. Hence, in addition to the SRIRMD sampling strategy, the HS [15], TEI [19], and EFDC [20] sampling strategies are also utilized to construct the surrogate models for the derivative functions of the state equation. The termination criterion working with those sampling strategies are the model accuracy at DTPs $\epsilon <0.001$ and the successive relative improvement of the objective function $\Delta J<0.1$, called MASRI. Hence, the BDCDO solving methods combined with those sampling strategies and MASRI are named HS-MASRI, TEI-MASRI, EFDC-MASRI, and SRIRMD-MASRI, respectively. To validate the robustness of these methods, each method is tested ten times. In each test, the initial samples for those methods are the same, the number of initial points $N_0=50$, and maximum number of new samples per iteration ${\Delta }_N=20$.

The test results of different methods are listed in

Table 3. The test results of the different methods in Example 1.

Method		HS-MASRI		TEI-MASRI		EFDC-MASRI		SRIRMD-MASRI		SRIRMD-STOR
Index		NoS	J	NoS	J	NoS	J	NoS	J	NoS	J
Test	1 2 3 4 5 6 7 8 9 10	350 350 350 350 350 350 350 350 350 350	5161.4 5186.3 5177.1 5166.0 5170.8 5178.0 5162.8 5184.7 5224.6 5258.2	346 240 325 319 315 345 346 267 346 346	5152.3 5157.5 5169.4 5156.4 5152.4 5154.3 5155.5 5152.6 5156.3 5153.3	185 249 230 225 238 215 225 231 197 230	5153.4 5159.7 5153.3 5155.0 5165.7 5155.0 5160.7 5157.6 5155.7 5161.9	153 207 193 183 177 184 150 152 145 150	5152.4 5153.0 5152.8 5152.6 5152.8 5152.7 5153.0 5152.3 5152.4 5152.5	130 170 150 159 159 167 130 130 121 130	5153.1 5153.6 5153.1 5152.6 5153.1 5153.0 5153.2 5151.8 5152.6 5153.0

, and the comparisons are depicted in Figure 4. NoS denotes the sample point size (or number of expensive valuations), and $J$ is the value of the performance index. From Table 1 and Figure 4, it can be observed that the HS-MASRI, TEI-MASRI, EFDC-MASRI, and SRIRMD-MASRI methods require an average of 350, 320, 223, and 170 samples to construct the surrogate models of the derivative functions, respectively. the average performance indexes are 5187.0, 5156.0, 5157.8, and 5152.7. While the SRIRMD-STOR method only needs 145 samples to construct the surrogate models, the average performance index is 5152.9. Compared to HS-MASRI, TEI-MASRI, and EFDC-MASRI, SRIRMD-STOR has the smallest error with the exact value. In SRIRMD-MASRI and SRIRMD-STOR, the sampling strategies are the same, and the termination criterion are different. SRIRMD-MASRI has a little higher solving accuracy than SRIRMD-STOR, while needing many more samples to construct the surrogate models, which indicates that the termination criterion for STOR is more effective than the termination criterion for MASRI by avoiding redundant iterations. As a result, for this numerical example, the BDCDO solving framework combined with SRIRMD and STOR can obtain a more robust and accurate approximate solution with less samples compared to the other methods.

To visualize the sampling outcomes of the HS-MASRI, TEI-MASRI, EFDC-MASRI, and SRIRMD-STOR methods, the distribution of sample points in different methods and the phase diagrams of optimal trajectories between some state variables are displayed in Figure 5 and Figure 6. The black dots are the initial sample points, the black stars are the new samples obtained via different sampling methods, and the red solid lines are the state trajectories optimized based on the surrogate models. From Figure 5 and Figure 6, it can be observed that the new sample points in the TEI, EFDC, and SRIRMD sampling strategies are located in the vicinity of the state trajectories, except for the HS strategy.

Figure 7 draws the iterative processes of the trajectories of partial state components in the SRIRMD-STOR method. As can be viewed from Figure 7, the trajectories of different state components converge gradually with the iterations.

Figure 8 records the convergence processes of the state component trajectory overlap ratio and the state trajectory overlap ratio in the SRIRMD-STOR method, ${\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6$, and $\mathrm{A}$ are the trajectory overlap ratios of the state components $S,T,L_1,L_2,I_1,I_2$, and state $\xi$, respectively. According to Figure 8, the state trajectory overlap ratio $\mathrm{A}$ converges to 1 with the convergences of all state components ${\alpha }_i$, and ${\alpha }_i\ge \mathrm{A}$, which verifies Theorem 1 and Theorem 2.

Moreover, to better compare the exact solution based on the mathematical model and the approximate solution based on the surrogate model, Figure 9 and Figure 10 show the comparison of the trajectories of the exact and approximate solutions about the state variables and control curves in Example 1, with the dashed lines being the exact solution and the thick solid lines being the approximate solution obtained by the SRIRMD-STOR method. Although the approximate and exact solutions of the control variables do not completely coincide, the trends remain consistent, and the effects on the trajectories of state variables are within acceptable limits.

4.2. Example 2: A Mathematical Nonlinear Dynamic Codesign and Optimization Problem

Example 2 is a dynamic co-design and optimization problem with one physical design parameter, one control input, and three state variables, the mathematical model of the problem is as follows

$$\begin{array}{ll}\underset{x_p,u(t)}{\min }& J=\sin (x_p)\cdot \cos (x_p)\cdot t_f\\ s.t.& \dot{x}=v\cdot \sin (u)\\ & \dot{y}=-v\cdot \cos (u)\\ & \dot{v}=10\cos (u)\\ & x_p\in [-\pi /2,\pi /2],u\in [-\pi /2,\pi /2],\\ & x\in [50,-50],y\in [0,-50],v\in [0,6]\end{array}$$

(24)

where only the objective function contains the plant design parameter $x_p$. The state variables $\xi$ include $\left[x,y,v\right]$, and the initial value and final value of $\xi$ are $\xi (t_0)=[0,0,0]$ and $\xi (t_f)=[2,-2,6]$, respectively. To solve this DOP based on the mathematical model, the maximum number of iterations of the NLP solver is set to 20, the solution accuracy is set to ${10}^{-6}$, the optimal plant design parameter is $-0.7854$, and the optimal performance index is $-50.00$. Similar to Example 1, the SRIRMD-STOR, TEI-MASRI, EFDC-MASRI, and SRIRMD-MASRI methods are utilized to solve this example.

The optimal physical design parameters and optimal objective values by different methods are listed in

Table 4. The computational cost, optimal plant parameters, and optimal objective values in Example 2.

	Mathematical Model	TEI-MASRI	EFDC-MASRI	SRIRMD-MASRI	SRIRMD-STOR
NoS	1255262	93	55	100	90
XP	−0.7854	−0.7854	1.5708	−0.7854	−0.7854
J	−50.00	−50.00	0.0000	−50.00	−50.00

. As can be observed from Table 2, the approaches based on the surrogate model significantly reduce the number of valuations of the mathematical model and save computational budget. The TEI-MASRI, SRIRMD-MASRI, and SRIRMD-STOR methods obtain the optimal physical parameters and objective values, and the SRIRMD-STOR method performs better in terms of efficiency with the assistant of the termination criterion STOR.

To intuitively demonstrate the sampling effect of the SRIRMD sampling strategy, the distribution of sample points and the phase diagram of the optimal trajectory between $y$ and $v$ are drawn in Figure 11. The left figure plots the distribution of all sample points (initial sample points and new sample points), and the right figure displays the positions of new samples, as well as a local zoomed-in view of the left figure near the optimal trajectory. It can be found from Figure 11 that the new sample points in the SRIRMD sampling strategy are all located in the vicinity of the state trajectory.

Figure 12 is the iteration processes of the trajectories of the state components $x$ and $v$ in the SRIRMD-STOR method. As it can be observed in Figure 12, the trajectories of different state components converge gradually with the iterations, and the trajectories of the 12th and 13th iterations tend to coincide. Meanwhile, Figure 13 exhibits the control curves obtained using the SRIRMD-based method.

Figure 14 records the convergence processes of the state component trajectory overlap ratio and the state trajectory overlap ratio in the SRIRMD-STOR method. ${\alpha }_1,{\alpha }_2,{\alpha }_3$ and $\mathrm{A}$ are the trajectory overlap ratios of the state components $x,y,v$, and state $\xi$, respectively. In light of Figure 14, the state trajectory overlap ratio $\mathrm{A}$ converges to 1 with the convergences of all ${\alpha }_i$, and ${\alpha }_i\ge \mathrm{A}$, which verifies Theorem 1 and Theorem 2.

5. Engineering Examples

In this section, the SRIRMD-STOR method is applied to optimize the 3-DOF Manutec r3 system and the horizontal axis wind turbine (HAWT) system. Different from the numerical examples, the derivative functions of the state equations are unknown in those dynamic systems, since only the simulation models can be accessed. In addition to the methods proposed in this work, other existing methods are also used to optimize the 3-DOF Manutec r3 system and HAWT system.

5.1. The BDCDO of 3-DOF Manutec r3 System

The industrial robot Manutec r3 [33], as shown in Figure 15, has six links, and only the three degrees of freedom associated with positioning are considered in this research, simplifying the robotic system to a three degrees of freedom (DOF) dynamic system. The goal of the 3-DOF robot co-design and optimization problem is to identify the optimal solution in the design space and state space so that the robot can move from the initial position to the specified position in the minimum time while satisfying the associated state equation constraints, upper and lower bound constraints of plant parameters, and state and control variables. The trajectory optimization formulation of the 3-DOF Manutec r3 system is described as follows:

$$\begin{array}{ll}\underset{x_p,u(t)}{\min }& J=t_f\\ s.t.& \dot{\xi }(t)=f(x_p,\xi (t),u(t),t)\\ & x_p\in [x_L,x_U]\\ & \xi (t)\in [{\xi }_L,{\xi }_U]\\ & u(t)\in [u_L,u_U]\end{array}$$

(25)

where the plant design parameters $x_p$ are the lengths of link 1 and link 2 $\left[L_1,L_2\right]$. The state variables $\xi$ consists of the relative angles of rotation $\left[\alpha ,\beta ,\gamma \right]$ and the relative angular velocities $[\dot{\alpha },\dot{\beta },\dot{\gamma }]$ between the connecting links, and the control variables $u$ includes the standardized torque controls $\left[u_1,u_2,u_3\right]$. The design intervals of the plant parameters $x_p$ are $L_1\in [0.4,0.5]$ and $L_2\in [0.9,1.0]$. The initial value and final value of $\xi$ are $\xi (t_0)=[-2,-2.5,-2,0,0,0]$ and $\xi (t_f)=[2,2.5,2,0,0,0]$, and the upper and lower bounds of $\xi$ are ${\xi }_U=[3,3,3,5,10,15]$ and ${\xi }_L=[-3,-3,-3,-5,-10,-15]$. The control variables are subject to interval constraints $u_U=[10,10,10]$ and $u_L=[-10,-10,-10]$. To solve this problem, the maximum number of iterations of the NLP solver is set to 20, and the solution accuracy is set to ${10}^{-6}$. In the original dynamic model, the physical design parameters are $x_p=[0.4500,0.9500]$, and the optimal solution yields the objective value as 0.9082. As with the above example, in addition to the SRIRMD-STOR method, the TEI-MASRI and EFDC-MASRI methods are also applied to optimize the 3-DOF Manutec r3 system. In the SRIRMD-STOR method, the number of initial points is $N_0=25$, and the maximum number of new samples per iteration is ${\Delta }_N=20$. It is notable that $\epsilon <0.001$ and $\Delta J<0.0001$ in the termination criterions of the TEI-MASRI and EFDC-MASRI methods.

The optimal physical design parameters and optimal objective values by different methods are presented in

Table 5. The computational cost, optimal plant parameters, and optimal objective values in the 3-DOF Manutec r3 system.

	Dynamic Model	TEI-MASRI	EFDC-MASRI	SRIRMD-MASRI	SRIRMD-STOR
NoS	3422	215	256	228	203
[L1,L2]	[0.4500, 0.9500]	[0.4151, 1.0000]	[0.4121, 0.9339]	[0.4019, 0.9982]	[0.4019,0.9982]
J	0.9082	0.9067	0.9064	0.9060	0.9060

. It can be observed from Table 5 that the surrogate-model-based approaches greatly reduce the number of valuations of the dynamic system and save computational costs compared to the original dynamic-model-based approach. Meanwhile, compared with the original physical design parameters, the optimal physical design parameters obtained by the TEI-MASRI, EFDC-MASRI, SRIRMD-MASRI, and SRIRMD-STOR methods shorten the working time to complete the specified task to a certain extent and improve the efficiency and performance of the robot arm. Among those surrogate model-based methods, the SRIRMD-STOR method obtains the best performance metrics using less samples. Therefore, the BDCDO solving framework combined with SRIRMD and STOR is a better alternative for the co-design and optimization problem of the 3-DOF Manutec r3 system.

To graphically demonstrate the sampling outcome of the SRIRMD sampling strategy, the distribution of sample points and the phase diagram of the optimal trajectory between different state variables are shown in Figure 16. The black dots are the initial samples, the black stars are the new samples obtained via SRIRMD, and the red solid lines are the state trajectories optimized based on the surrogate models. As Figure 16 reveals, the new sample points in the SRIRMD sampling strategy are all situated nearby the state trajectories.

Figure 17 graphs the trajectory iteration processes for the state components $\alpha ,\beta ,\gamma ,\dot{\alpha },\dot{\beta }$, and $\dot{\gamma }$ in the SRIRMD-STOR method. As visible in Figure 17, the trajectories of the different state components converge gradually as the iterations proceed, and the trajectories of the 17th and 18th iterations tend to coincide. Meanwhile, Figure 18 exhibits the control curves obtained by the SRIRMD method.

Figure 19 records the convergence processes of the state component trajectory overlap ratio and the state trajectory overlap ratio in the SRIRMD-STOR method, ${\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5,{\alpha }_6$, and $\mathrm{A}$ are the trajectory overlap ratios of the state components $\alpha ,\beta ,\gamma ,\dot{\alpha },\dot{\beta },\dot{\gamma }$, and state $\xi$, respectively. According to Figure 19, the state trajectory overlap ratio $\mathrm{A}$ converges to 1 with the convergences of all ${\alpha }_i$, and ${\alpha }_i\ge \mathrm{A}$, which verifies Theorem 1 and Theorem 2.

5.2. The BDCDO of the Horizontal Axis Wind Turbine (HAWT)

The design optimization problem of the HAWT system [4,15,20] is a complex co-design BDOP involving structural parameters and control variables, which can be simulated and estimated by the Advanced Wind Turbine program blade 27 (AWT27) in the Open FAST project. As shown in Figure 20, the structural parameters $x_p$ in this paper mainly include the hub radius $R_h$, blade length $L_b$, and tower height $H_t$. The control variable $u$ is the generator torque $G_t$. The fore-aft tower-top displacement ${\xi }_1$, the side-to-side tower-top displacement ${\xi }_2$, the fore-aft tower-top velocity $v_1$, the side-to-side tower-top velocity $v_2$, and rotor speed $\omega$ are regard as the state variables $\xi$. The co-design formulation of the HAWT system is expressed as follows:

$$\begin{array}{ll}\underset{x_p,u(t)}{\min }& J=w_1{m}_s(x_p)+w_2{\displaystyle {\int }_0^{t_f}{(\lambda (t)-{\lambda }_{*}(t))}^{2}dt}\\ s.t.& \dot{\xi }(t)=f(\xi (t),x_p,u(t),t)\\ & {\Vert {\xi }_1(t)\Vert }_{\infty }-{\xi }_{1\max }\le 0\\ & {\Vert {\xi }_2(t)\Vert }_{\infty }-{\xi }_{2\max }\le 0\\ & P(V_e)-P_{e\min }\ge 0\\ & x_p\in \left[x_L,x_U\right]\end{array}$$

(26)

where $m_s(x_p)$ is the mass of the wind turbine, ${\displaystyle {\int }_0^{t_f}{(\lambda (t)-{\lambda }_{*}(t))}^{2}dt}$ is the sum of the deviations of the wind blade tip tangential velocity ratio $\lambda (t)$ and the optimal velocity ratio ${\lambda }_{*}(t)$ over a period of time, $\lambda$ is the ratio of leaf tip tangential velocity $\omega \cdot L_b$ to wind speed $v$, ${\lambda }_{*}$ can be calculated by the power coefficient function, and $w_1$ and $w_2$ are the weights of the two terms, respectively. Therefore, the optimization objective of the system is to minimize the sum of the mass and the deviations of speed ratios. What is more, the HAWT system needs to satisfy the structural deflection constraints and ${\Vert {\xi }_2(t)\Vert }_{\infty }$. When the wind speed reaches the rated wind speed $V_e$, the system has to reach the minimum rated power $P_{e\min }$. Since the simulation of the HAWT system involves several disciplines, the RHS function $\dot{\xi }(t)=f(\xi (t),X_p,u(t),t)$ of the system is highly nonlinear, and AWT27 takes several seconds to execute a simulation valuation.

The co-design and optimization problem of HAWT is optimized based on the finite difference technique without using the surrogate model, and the standard optimal objective 805.4801, hub radius $R_h=1.2000$, blade length $L_b=13.7330$, and tower height $H_t=32.3944$ are obtained. Qiao et al. [20] adopted the HS-MASRI, EFDC-MASRI, and TEI-MASRI methods to solve this problem, while the SRIRMD method is used to build the surrogate model of a derivative function in the HAWT system. In the SRIRMD-STOR method, the number of initial points $N_0=100$, and maximum number of new samples per iteration ${\Delta }_N=10$. The optimization outcomes of those different methods are listed in

Table 6. The computational cost, optimal plant parameters, and optimal objective values in the HAWT system.

	Original System	HS-MASRI	EFDC-MASRI	TEI-MASRI	SRIRMD-MASRI	SRIRMD-STOR
NoS	546400	540	291	460	240	210
[Rh Lb Ht]	[1.2000,13.7330,32.3944]
J	805.4801	806.0783	806.3677	807.1469	805.3558	805.3550
Absolute Error	0	0.5982	0.8876	1.6668	0.1243	0.1251

. Obviously, the original system-based optimization solution method is costly and necessitates extensive simulation evaluations of AWT27. In contrast, the surrogate model-based optimization methods significantly reduce the number of running valuations of AWT27 and decrease the computational costs. In these model-based optimization solutions, the plant parameters converge to the standard solution $[1.2000,13.7330,32.3944]$. Due to the accuracy of the surrogate models, the objective values obtained by various methods are different. Nevertheless, the errors with the standard solution are within the allowed range. More importantly, it is clear that the SRIRMD-STOR method uses the least number of samples in addressing the co-design problem of the HAWT system while it has the higher solution accuracy, improving the solution efficiency and conserving the computational resources.

The wind speed curve input to the HAWT system for a certain time period is displayed in Figure 21. The evolution of the generator torque $G_t$, the fore-aft tower-top displacement ${\xi }_1$, the side-to-side tower-top displacement ${\xi }_2$, and the rotor speed $\omega$ with this wind speed curve are shown in Figure 22. As can be observed from the figure, in order to optimize the HAWT system, the trend of the control and state variables obtained from the optimization solution is highly consistent with the trend of the wind speed.

Figure 23 records the convergence processes of the state component trajectory overlap ratio and the state trajectory overlap ratio in the SRIRMD-STOR method, and ${\alpha }_1,{\alpha }_2,{\alpha }_3,{\alpha }_4,{\alpha }_5$, and $\mathrm{A}$ are the trajectory overlap ratios of the state components ${\xi }_1,{\xi }_2,v_1,v_2,\omega$, and state $\xi$, respectively. According to Figure 23, the state trajectory overlap ratio $\mathrm{A}$ converges to 1 with the convergences of all ${\alpha }_i$ in the $11th$ iteration. ${\alpha }_i\ge \mathrm{A}$, which verifies Theorem 1 and Theorem 2.

6. Conclusions

The BDCDO solving framework based on the surrogate models of the derivative functions in the state equation is an effective approach to solve the black-box dynamic co-design and optimization problem. For efficient construction of the surrogate models for the right hand-side functions of the state equation, a novel adaptive sequential sampling strategy, called SRIRMD, was proposed in this work. This strategy refines the surrogate models by selecting suitable sample points from the trajectory discrete points. At the same time, to quantify the convergence and intuitively reflect the convergence trend of the solution during the solving process, a new termination criterion, called the state trajectory overlap ratio (STOR), was also introduced. Finally, the BDCDO solving framework combined with SRIRMD and STOR was utilized to address two numerical optimization problems and two engineering co-design and optimization problems. The numerical examples indicate that the SRIRMD sampling strategy proposed in this work was superior to the existing sampling strategies with respect to both the solution accuracy and robustness. The 3-DOF robot co-design and optimization problem and the horizontal axis wind turbine co-design and optimization problem revealed that the BDCDO solving framework combined with SRIRMD and STOR is a feasible and efficient tool to optimize the black-box dynamic systems. In summary, the proposed sampling strategy and the termination criterion in this research not only improve the efficiency of the BDCDO solving framework but also save the computational budget.