**The SRIRMD-STOR Method**

*Step 1:* Apply LHS method to sample initial points set *S*<sup>0</sup> in the design domain consisting of feasible regions of physical design parameters **x***p*, state variables ξ and control inputs **u**.

*Step 2:* Construct the initial surrogate model ˆ *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** (0) *<sup>p</sup>* , **Ξ**(0), **Θ**(0), ˆ *f* (0)] and obtain the current optimal plant design

parameters, state trajectories, control curves, and performance index [**x** (1) *<sup>p</sup>* , **Ξ**(1), **Θ**(1), *J*(1)]. *Step 4:* Calculate the state component trajectory overlap ratios *α<sup>i</sup>* of all state variables according to the initial guess trajectories and the current optimal trajectories, then calculate the state trajectory overlap ratio A. If A > A0, terminate the solving process; otherwise, go to Step 5. *Step 5:* Employ the SRIRMD strategy to select new samples *xnew* from the current DTPs, update

the samples set *S*1, and rebuild the surrogate model ˆ *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** (*l*) *<sup>p</sup>* , **Ξ**(*l*), **Θ**(*l*), ˆ *f* (*l*)]; and acquire the latest optimal plant parameters, state trajectories, control curves, and performance index. Note: *l* starts from 1. *Step 8:* Calculate *α<sup>i</sup>* and A. If A > A0, stop the solving process; otherwise, go to Step 5.
