[insert Figure 3 here]
The computational difficulty of the MINLP problem depends on the number of ingredient candidates and the complexity of the adopted models. A large number of ingredient candidates often create a combinational problem. Many rigorous models (e.g., thermodynamic and transport phenomena models) involve nonlinear and nonconvex equations. Along with complex surrogate models (e.g., neural network), the optimization problem is prone to convergence failure if the problem is directly solved using standard MINLP solvers. Some of these problems can be handled with advanced algorithms. For instance, Schweidtmann and Mitsos42,43 recently developed and applied an efficient global solver for ANN embedded MINLP problems. Alternatively, the problem can be resolved by reformulating the optimization problem.44 Two techniques are proposed for enhancing optimization convergence and finding better solutions (Figure 3). If the problem can be directly solved, the solution is sent for experimental validation. Otherwise, generalized disjunctive programming (GDP) can be used because of the need to calculate the multiple intermediate variables in the mechanistic models. If the GDP problem still cannot be solved or better solutions are needed, the model(s) in use are replaced with alternative model(s) and repeat the calculations.

GDP reformulation

As can be seen in Eq. 11, even if i -th ingredient candidate is not selected, its intermediate variables (\(IM_{i}^{m}\)) must be calculated. Also, forcing \(V_{i}\) to 0 may lead to singularity at\(V_{i}=0\) for some models (e.g., logarithmic function). Thus, when mechanistic models are employed and multiple intermediate variables are calculated through complex equations, the redundant constraints and singularities can lead to convergence failure. Similar to the tray selection problem in distillation column design,44cosmetic formulation problem can be formulated using GDP. As an alternative way to program discrete-continuous problem, GDP is a logic-based method containing Boolean and continuous variables. Constraints are expressed as disjunctions, algebraic equations, and logic propositions.44 The following disjunction can be used to express the bounds and intermediate variables for Eq. 11.
\(\par \begin{bmatrix}Y_{i}\\ VL_{i}\leq V_{i}\leq VU_{i}\\ \par \begin{matrix}IM_{i}^{m}=IMG^{m}(V_{i},ms)\\ \end{matrix}\\ \end{bmatrix}\bigvee\par \begin{bmatrix}\neg Y_{i}\\ V_{i}=0\\ \par \begin{matrix}IM_{i}^{m}=0\\ \end{matrix}\\ \end{bmatrix}\),\(i\in\left\{I_{A,1},I_{A,2},\ldots,I_{Z,z}\right\}\) (12)
where \(Y_{i}\) is the Boolean variable for ingredient selection. If thei -th ingredient is selected, the bounds on \(V_{i}\) are fulfilled and its intermediate variables \(IM_{i}^{m}\) are calculated. Otherwise, they are not calculated and simply set as 0.
To solve a GDP problem, it is often transformed back into MINLP using big-M or convex-hull relaxation to take advantage of standard MINLP solvers. It is found that the big-M method is more appropriate in solving mixture design problem since singularity issue can still occur in the convex-hull relaxation.12 After transforming the above disjunction using big-M approach, the cosmetic formulation problem is reformulated below. \(S_{i}\) has a one-to-one correspondence with \(Y_{i}\). bm is a sufficiently large parameter.
\(\operatorname{}{q=f(V_{I_{A,1}},V_{I_{A,2}},\ldots,V_{I_{Z,z}},ms)}\)(13)
s.t. \(PL^{k}\leq P^{k}\leq PU^{k}\), \(k\in K=MM\cup SM\)
\(\left\{\par \begin{matrix}VL_{i}-bm\bullet(1-S_{i})\leq V_{i}\leq VU_{i}+bm\bullet(1-S_{i})\\ -bm\bullet S_{i}\leq V_{i}\leq bm\bullet S_{i}\\ \par \begin{matrix}\text{IM}G^{m}\left(V_{i},ms\right)-bm\left(1-S_{i}\right)\leq IM_{i}^{m}\leq IMG^{m}\left(V_{i},ms\right)+bm(1-S_{i})\\ -bm\bullet S_{i}\leq IM_{i}^{m}\leq bm\bullet S_{i}\\ \end{matrix}\\ \end{matrix}\right.\ \) Big-M constraint
\(P^{m}=G^{m}(IM_{I_{A,1}}^{m},IM_{I_{A,2}}^{m},\ldots,IM_{I_{Z,z}}^{m})\),\(m\in MM\)
\(P^{s}=g^{s}(V_{I_{A,1}},V_{I_{A,2}},\ldots,V_{I_{Z,z}},ms)\),\(s\in SM\)
\(H\left(S_{I_{A,1}},S_{I_{A,2}},\ldots,S_{I_{Z,z}},V_{I_{A,1}},V_{I_{A,2}},\ldots,V_{I_{Z,z}}\right)\leq 0\)
\(msL\leq ms\leq msU\), \(S_{i}\in\left\{0,1\right\}^{i}\),\(\sum_{i}{V_{i}=1}\),\(i\in\left\{I_{A,1},I_{A,2},\ldots,I_{Z,z}\right\}\)

Model substitution

Some rigorous mechanistic models are too complicated to be directly used for optimization even if they are programmed using GDP. There is always a trade-off between model accuracy and traceability. In this case, the complicated but accurate rigorous models can be replaced by simple short-cut model or surrogate model to reduce the computational effort and to seek out even better solutions. A surrogate model is a good choice when it is relatively easy to generate simulation data as training data from the rigorous model. Although the model accuracy is reduced, it is easier to solve and to obtain the global solution.42,45–47 After model substitution, the newly generated optimization problem should be solved and the optimal solution obtained can be denoted as \(V_{i}^{*}\). This solution must be validated using the original rigorous mechanistic models. If the validation fails, the newly generated optimization problem should be re-solved by adding the equation below to remove this solution (\(V_{i}^{*}\)) that fails validation. Otherwise, the solution can be sent for experimental verification.
\(\sum_{i}\left(V_{i}^{**}-V_{i}^{*}\right)^{2}\geq tol\) (14)
Here, \(V_{i}^{**}\) is the solution of a new round optimization. The parameter tol is a small tolerance.