*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 & f = \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\\ \text{s.t.} & \dot{S} = \Lambda - \frac{13}{30000}SI\_{1} - \frac{0.029}{50000}SI\_{2} - 0.0143S\\ & \dot{T} = 2u\_{1}L\_{1} - 0.0143T + (1 - 0.5(1 - u\_{2}))I\_{1} - \frac{13}{30000}TI\_{1} - \frac{0.029}{50000}TI\_{2} \\ & \dot{L}\_{1} = \frac{13}{30000}SI\_{1} - 0.6143L\_{1} - 2u\_{1}L\_{1} + 0.4(1 - u\_{2})I\_{1} + \frac{13}{30000}TI\_{1} - \frac{0.029}{50000}L\_{1}I\_{2} \\ & \dot{L}\_{2} = 0.1(1 - u\_{2})I\_{1} - 0.0143L\_{2} + \frac{0.029}{50000}(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} + L\_{2} + T = 30000 \end{array} \tag{23}$$

where the state variables ξ = [*S*, *T*, *L*1, *L*2, *I*1, *I*2] *<sup>T</sup>*. The upper and lower bounds of control inputs **u** = [*u*1, *u*2] *<sup>T</sup>* are **<sup>u</sup>***<sup>U</sup>* = [0.95, 0.95] *<sup>T</sup>* and **<sup>u</sup>***<sup>L</sup>* = [0.05, 0.05] *<sup>T</sup>*. The finial time *tf* 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−<sup>4</sup> 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 *ε* < 0.001 and the successive relative improvement of the objective function Δ*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*<sup>0</sup> = 50, and maximum number of new samples per iteration Δ*<sup>N</sup>* = 20.

The test results of different methods are listed in Table 3, 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.


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

**Figure 4.** Box-plots of test results of different methods in Example 1.

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 Figures 5 and 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 Figures 5 and 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, *α*1, *α*2, *α*3, *α*4, *α*5, *α*6, and A are the trajectory overlap ratios of the state components *S*, *T*, *L*1, *L*2, *I*1, *I*2, and state ξ, respectively. According to Figure 8, the state trajectory overlap ratio A converges to 1 with the convergences of all state components *αi*, and *α<sup>i</sup>* ≥ A, which verifies Theorem 1 and Theorem 2.

**Figure 5.** The phase diagrams of optimal trajectories and distribution of samples between state variables *S* and *T*.

**Figure 6.** The phase diagrams of optimal trajectories and distribution of samples between state variables *I*<sup>1</sup> and *I*2.

**Figure 7.** Trajectory iterative processes of state components *S*, *L*1, *L*2, and *I*1.

**Figure 8.** The convergence processes of A and all *α<sup>i</sup>* in the numerical example 1.

Moreover, to better compare the exact solution based on the mathematical model and the approximate solution based on the surrogate model, Figures 9 and 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.

**Figure 9.** The exact solution and approximate solution of the state variables.

**Figure 10.** The exact solution and approximate solution of the control variables.

## *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}\min\_{\mathbf{x}\_p, \mathbf{u}(t)} & J = \sin(\mathbf{x}\_p) \cdot \cos(\mathbf{x}\_p) \cdot t\_f \\ \text{s.t.} & \dot{\mathbf{x}} = \boldsymbol{v} \cdot \sin(\mathbf{u}) \\ & \dot{y} = -\boldsymbol{v} \cdot \cos(\mathbf{u}) \\ & \dot{v} = 10 \cos(\mathbf{u}) \\ & x\_p \in [-\pi/2, \pi/2], \boldsymbol{u} \in [-\pi/2, \pi/2], \\ & x \in [50, -50], y \in [0, -50], v \in [0, 6] \end{array} \tag{24}$$

where only the objective function contains the plant design parameter **x***p*. The state variables ξ include [*x*, *y*, *v*], and the initial value and final value of ξ are ξ(*t*0)=[0, 0, 0] and ξ(*tf*)=[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. 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.


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

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. *α*1, *α*2, *α*<sup>3</sup> and A are the trajectory overlap ratios of the state components *x*, *y*, *v*, and state ξ, respectively. In light of Figure 14, the state trajectory overlap ratio A converges to 1 with the convergences of all *αi*, and *α<sup>i</sup>* ≥ A, which verifies Theorem 1 and Theorem 2.

**Figure 11.** The phase diagrams of the optimal trajectory and distribution of the samples between state variables *y* and *v*.

**Figure 12.** Trajectory iterative processes of state components *x* and *v*.

**Figure 13.** The solution of the control variable *u*.

**Figure 14.** The convergence processes of A and all *α<sup>i</sup>* in the numerical example 2.
