**4. Proposed Solution Method**

In this paper, the interior point method (IPM) [29] is conducted to solve (25). Towards this end, the *Jacobian* and *Hessian* matrix of the trained DBN model is analyzed.

*4.1. Interior Point Method*

For the sake of simplification, model Equation (23) is firstly reformed as the following canonical form:

$$\begin{array}{c} \text{Max } F(\mathbf{x}), \\ \text{s.t. } G(\mathbf{x}) = \mathbf{0}, \\ \left[ H(\mathbf{x}) \le \mathbf{0} \right] = \left[ H\_{\boldsymbol{\epsilon}}(\mathbf{x}) \le \mathbf{0}, H\_{\boldsymbol{S}}(\mathbf{x}) \le \mathbf{0} \right] \end{array} \tag{24}$$

where *F*(*x*) is the objective function; *G*(*x*)=[*G*1(*x*), ... ,*G*m(*x*)]*<sup>T</sup>* is the nonlinear equality constraints; *H*(*x*)=[*H*1(*x*), ... ,*H*r(*x*)]*<sup>T</sup>* is the non-linear inequality constraints, and *Hc*(*x*), *HS*(*x*) are the constraints in (4) and the surrogate model in (23b), respectively; *r*,*m* are the number of inequality and equality constraints.

Use IPM to solve Equation (24), and the steps are as follows [29]:


$$\begin{aligned} L &= F(\mathbf{x}) - \boldsymbol{\zeta}^T \mathbf{G}(\mathbf{x}) - \boldsymbol{z}^T [H(\mathbf{x}) - \boldsymbol{l} - H\_{\text{min}}] - \boldsymbol{\omega}^T [H(\mathbf{x}) + \boldsymbol{u} - H\_{\text{max}}] \\ &- \mu \boldsymbol{\Sigma}^{\tau}\_{\,i=1} \log(l\_i) - \mu \boldsymbol{\Sigma}^{\tau}\_{\,i=1} \log(l\_i), \end{aligned} \tag{25}$$

where ζ,*z*,*ω*, are Lagrangian multipliers, respectively.

4. Calculate *μ* via Equation (26):

$$
\mu = \sigma (\mathbf{l}^T \mathbf{z} - \mathbf{u}^T \boldsymbol{\omega}) / 2\mathbf{r},\tag{26}
$$

where *σ* denotes the central parameter.

5. Consider the Karush–Kuhn–Tucker (KKT) conditions and adopt the Newton method, the matrix form of the modified equations can be deduced as:

$$\begin{array}{c} \mathbf{A} \cdot \Delta \mathbf{x} + (\partial G(\mathbf{x})/\partial \mathbf{x}) \cdot \Delta \mathbf{x} = \Phi, \\ (\partial G(\mathbf{x})^T/\partial \mathbf{x}) \cdot \Delta \mathbf{x} = \mathbf{G}, \end{array}$$

$$\begin{array}{c} \mathbf{A} = (\partial^2 G(\mathbf{x})/\partial \mathbf{x}^2)\mathbf{\hat{\zeta}} + (\partial^2 H(\mathbf{x})/\partial \mathbf{x}^2)(\mathbf{z} + \boldsymbol{\omega}) - (\partial^2 F(\mathbf{x})/\partial \mathbf{x}^2) \\ + (\partial H(\mathbf{x})/\partial \mathbf{x})(\mathbf{u}^{-1}\boldsymbol{\omega} - \boldsymbol{l}^{-1}\mathbf{z})(\partial H(\mathbf{x})/\partial \mathbf{x})^T, \\ \boldsymbol{\Phi} = -L\_{\mathbf{x}} - (\partial H(\mathbf{x})/\partial \mathbf{x})[L^{-1}(L\_{\mathbf{l}}\boldsymbol{\mu} + \mathbf{Z}L\_{\mathbf{z}}) + \mathbf{U}^{-1}(L\_{\mathbf{u}}\boldsymbol{\mu} + \mathbf{W}L\_{\mathbf{u}})] \end{array} \tag{27}$$

where *Lx*, *Lz*, *Lω*, *Ll <sup>μ</sup>* and *L<sup>u</sup> <sup>μ</sup>* are the partial derivatives of *L* to *x*, *z*, *ω*, *l* and *u*. Besides, the corrections of *z*, *ω*, *l* and *u* can be calculated via (28):

$$\begin{cases} \Delta \mathbf{z} = \mathbf{L}^{-1} L\_l \mu - \mathbf{L}^{-1} \mathbf{Z} \Delta l, \\ \Delta l = (\partial H(\mathbf{x}) / \partial \mathbf{x})^T \Delta \mathbf{x} - L\_{z\prime} \\ \Delta \boldsymbol{\omega} = \mathbf{U}^{-1} L\_{\mu} \mu - \mathbf{U}^{-1} \mathbf{W} \Delta \boldsymbol{\omega}, \\ \Delta \boldsymbol{\mu} = - (\partial H(\mathbf{x}) / \partial \mathbf{x})^T \Delta \mathbf{x} + L\_{\omega} \end{cases} \tag{28}$$

where *Z* = diag(*z*); *W* = diag(*ω*); *L* = diag(*l*); and *U* = diag(*u*).

6. Use the corrections calculated via Equations (27) and (28) to update the variables as follows:

$$\begin{array}{l} \mathbf{x}^{(k+1)} = \mathbf{x}^{(k)} + \alpha\_p \Delta \mathbf{x}, \; \mathsf{L}^{(k+1)} = \mathsf{L}^{(k)} + \alpha\_d \Delta \mathsf{L}\_r \\\ I^{(k+1)} = I^{(k)} + \alpha\_p \Delta I, \; \mathsf{z}^{(k+1)} = \mathsf{z}^{(k)} + \alpha\_d \Delta \mathsf{z}\_r \\\ \mathsf{u}^{(k+1)} = \mathsf{u}^{(k)} + \alpha\_p \Delta \mathsf{u}, \; \mathsf{w}^{(k+1)} = \mathsf{w}^{(k)} + \alpha\_d \Delta \omega\_r \end{array} \tag{29}$$

where the step size α*<sup>p</sup>* and α*<sup>d</sup>* for each update are shown in Equation (30).

$$\begin{aligned} \alpha\_p &= 0.9995 \text{min[min(-I\_i/\Delta I\_i, I\_i < 0; -I\_i/\Delta \mu\_i, \mu\_i < 0), 1]}, \\ \alpha\_d &= 0.9995 \text{min[min(-z\_i/\Delta z\_i, z\_i < 0; -\omega\_i/\Delta \omega\_i, \omega\_i > 0), 1], i = 1, \dots, r \end{aligned} \tag{30}$$

7. Termination condition: if (*l <sup>T</sup><sup>z</sup>* <sup>−</sup> *<sup>u</sup>Tω*) < *<sup>ε</sup>*, the current *<sup>x</sup>* is output; else re-execute (4) to (6). Here *ε* represents the specified threshold.
