[insert Figure 4 here]

LD50

\(LD_{50}\) of the perfume solution is calculated by Eq. 24. It depends on the toxicity of ingredients (\(LD_{50,i}\)) and the mass fraction (\(m_{i}\)) converted from the volume fraction \(V_{i}\).
\(LD_{50}=\frac{1}{\sum_{i=1}^{48}\frac{m_{i}}{LD_{50,i}}}\) (24)
\(m_{i}=\frac{V_{i}\bullet\rho_{i}}{\sum_{j=1}^{48}{V_{j}\bullet\rho_{j}}}\)(25)
where \(LD_{50,i}\) and the density (\(\rho_{i}\)) for the 48 ingredient candidates are given in Table S2.

Flash point

The flash point (\(T_{\text{fp}}\)) of a flammable liquid mixture can be theoretically determined based on the Le Chatelier’s mixing rule.53
\(\sum_{i=1}^{48}\frac{\text{FP}P_{i}}{\text{FPLF}L_{i}}=1\) (26)
where \(\text{FP}P_{i}\) and \(\text{FPL}FL_{i}\) are the partial pressure and lower flammability limit of the i -th ingredient candidate at the flash point, respectively. \(\text{FPL}FL_{i}\) is calculated by
\(\text{FPL}FL_{i}=LFL_{i}^{*}-\frac{0.182\times(T_{\text{fp}}-298)}{Hc_{i}}\)(27)
where \(Hc_{i}\) and \(\text{LF}L_{i}^{*}\) are the heat of combustion and lower flammability limit at 298 K (see Table S2), respectively.\(\text{FP}P_{i}\) is calculated via the vapor-liquid equilibrium in Eq. 28. The UNIFAC model is used to calculate the activity coefficient\(\text{FPγ}_{i}\) at the flash point. The mole fraction \(x_{i}\) is converted from mass fraction \(m_{i}\). \(\text{FPPsa}t_{i}\) is the saturated vapor pressure at flash point, which is calculated using the Antoine equation in Eq. 31.
\(\text{FP}P_{i}\ =\text{FPγ}_{i}\bullet x_{i}\bullet FPP\text{sat}_{i}\)(28)
\(\text{FPγ}_{i}=f_{\text{unifac}}\left(x_{i},T_{\text{fp}}\right)\)(29)
\(x_{i}=\frac{m_{i}}{MW_{i}\bullet\sum_{j}\frac{m_{j}}{MW_{j}}}\)(30)
\(\operatorname{}{\text{FPPsa}t_{i}}=A_{i}-\frac{B_{i}}{C_{i}+T_{\text{fp}}}\)(31)
The molecular weight \(MW_{i}\), UNIFAC parameters, and Antoine coefficients \(A_{i}\), \(B_{i}\), and \(C_{i}\) for the 48 ingredient candidates are given in Table S2.

Homogeneous solution

To ensure a homogeneous solution, the volume of selected organic fragrances must be less than their volume solubility (\(SV_{i,ew}\)) in the ethanol-water solvent system.
\(\frac{V_{i}}{(V_{47}+V_{48})}\leq SV_{i,ew},\ \ \ i=1,\ldots,46\)(32)
It is found that it is quite hard to calculate \(SV_{i,ew}\) using rigorous thermodynamic models due to the many missing parameters. In the literature, several short-cut models have been developed to predict\(SV_{i,ew}\). The log-linear mixture rule below is widely used.54
\(\log{\text{SV}_{i,ew}=}\log{SV_{i,w}}+\beta\bullet\log\frac{SV_{i,e}}{SV_{i,w}},\ \ \ \ i=1,\ldots,46\)(33)
\(\beta=\frac{V_{47}}{V_{47}+V_{48}}\) (34)
\(\log\frac{SV_{i,e}}{SV_{i,w}}=M\bullet\log K_{ow,i}+N,\ \ \ \ i=1,\ldots,46\)(35)
where \(SV_{i,e}\) and \(SV_{i,w}\) are the volume solubility in ethanol and water, respectively.\(\ K_{ow,i}\) is the n-octanol/water partition coefficient of the i -th candidate. \(M\) and \(N\) are the cosolvent constants. Based on experimental data, their values have been regressed as 0.81 and 0.85, respectively.

Odor type in top note

The fragrance molecules in a perfume solution first evaporate into the air through the liquid-gas interface. Then, the molecules diffuse in the air (assumed to be stagnant) and are detected at certain distance away. The processes of evaporation, diffusion, and detection have been modelled using chemical engineering principles and psychophysics.38,52,55 Perfume evaporation is simulated using Eq. 36 with an initial condition. The liquid molar changes are equal to the moles of ingredients transported through the interface (i.e., \(z=0\)).
\(\frac{dn_{i,t}}{\text{dt}}=C_{T}\bullet D_{i}\bullet A_{\lg}\ \bullet\left.\ \frac{\partial y_{i,t,z}}{\partial z}\right|_{z=0}\)(36)
Initial condition: \(n_{i,t=0}=n_{p}\bullet x_{i}\)
After discretization, Eq. 37 is obtained.
\(\frac{n_{i,t+t}-n_{i,t}}{t}=C_{T}\bullet D_{i}\bullet A_{\lg}\ \bullet\frac{{y_{i,t,z=z_{1}}-y}_{i,t,z=0}}{z_{1}}\)(37)
where \(n_{p}\) is the initial number of moles of perfume solution.\(C_{T}=P/RT\) is a constant.\(\ D_{i}\) and \(A_{\lg}\) are the diffusivity of i -th candidate and interfacial area, respectively.\(t\) and \(z_{1}\)are the time interval and the first distance interval, respectively. These parameters are given in Table S2.\(n_{i,t}\) is the number of moles of the i -th candidate in the liquid at time \(t\). \(y_{i,t,z}\) is the molar fraction of i -th ingredient candidate in the air at time \(t\) at distance z . It is calculated via vapor-liquid equilibrium.
\(y_{i,t,z=0}=\gamma_{i,t}\bullet x_{i,t}\bullet\frac{\text{Psa}t_{i}}{P}\)(38)
\(\gamma_{i,t}=f_{\text{unifac}}(x_{i,t},T_{r})\) (39)
\(x_{i,t}=\frac{n_{i,t}}{\sum_{i=1}^{48}n_{i,t}}\) (40)
where \(\gamma_{i,t}\) and \(x_{i,t}\) are the activity coefficient and mole fraction of i -th ingredient candidate at time t , respectively. \(\text{Psa}t_{i}\) is the saturated vapor pressure at room temperature \(T_{r}=298\ K\).
After evaporation, fragrance diffusion is modelled based on Fick’s 2nd law of diffusion with one initial condition and two boundary conditions (Eq. 41).
\(\frac{\partial y_{i,t,z}}{\partial t}=D_{i}\bullet\frac{\partial^{2}y_{i,t,z}}{\partial z^{2}}\)(41)
Initial condition: \(y_{i,t=0,z}=0\)
Boundary conditions: Eq. 38, \(y_{i,t,z=z_{\max}}=0\)
The initial condition assumes that no fragrances exist in the air before diffusion begins (i.e., \(t=0\)). The boundary conditions indicate that vapor-liquid equilibrium is maintained at the interface at any time (i.e., Eq. 38) and no fragrances exist beyond the maximum distance (\(z_{\max}=2m\)). This model is discretized using a non-uniform distance grid (Table S2) for reducing the computational difficulty. After discretization, we get
\(\frac{y_{i,t+t,z}-y_{i,t,z}}{t}=D_{i}\bullet\frac{\frac{y_{i,t,z+z_{j+1}}-y_{i,t,z}}{z_{j+1}}-\frac{y_{i,t,z}-y_{i,t,z-z_{j}}}{z_{j}}}{0.5\times(z_{j+1}+z_{j})}\),\(z\in[0,z_{\max}]\) (42)
where \(z_{j}\) and \(z_{j+1}\) are the distance intervals, respectively.
Any fragrance with a different concentration leads to a different intensity. Many theoretical models (e.g., Weber-Fenchner law, power law, and linear law) have been proposed for quantifying odor intensity. The power law is chosen here because it fits experimental data well. The intensity of the i -th odorant is defined as the ratio of its concentration in the air (\(c_{i}\) in g/m3) to its odor recognition threshold value (\(\text{OR}T_{i}\)), raised to a power\(oe_{i}\).52 With this, the odor intensity in the top note is determined based on the mole fraction of fragrances in the air at 5 minutes (\(t_{\text{tn}}\)) after application at a distance of 0.2 m (\(z_{\text{tn}}\)).
\(\psi_{i}=\left(\frac{c_{i}}{\text{OR}T_{i}}\right)^{oe_{i}}\) (43)
\(c_{i}=y_{i,t_{\text{tn}},z_{\text{tn}}}\bullet MW_{i}\bullet C_{T}\)(44)
Given multiple odorants, the one with the highest intensity is more strongly sensed and can be regarded as the major odor type. Thus, the dominant odor type in top note is expressed as
\(OTTN=i,\ \ if\ \psi_{i}=\psi_{\max}\operatorname{=}\left\{\psi_{i}\right\}\)(45)

Heuristics

Following Table 3, constraints for the Eau de parfum formulation are derived from dozens of modern Eau de parfum recipes.51It is found that the suggested number of ingredients for each fragrance note can be represented by Eq. 46-48. Eq. 49 shows that Eau de parfum usually contains 10-20% organic fragrances. The suggested volumetric proportions for top note and middle note are 15-25% and 30-40%, respectively (Eq. 50-51). The suggested volume fraction of water is 9-13%.49,52
\(3\leq\sum_{i=1}^{17}S_{i}\leq 6\) (46)
\(3\leq\sum_{i=18}^{33}S_{i}\leq 6\) (47)
\(2\leq\sum_{i=34}^{46}S_{i}\leq 5\) (48)
\(0.1\leq\sum_{i=1}^{46}V_{i}\leq 0.2\) (49)
\(0.15\bullet\sum_{i=1}^{46}V_{i}\leq\sum_{i=1}^{17}V_{i}\leq 0.25\bullet\sum_{i=1}^{46}V_{i}\)(50)
\(0.3\bullet\sum_{i=1}^{46}V_{i}\leq\sum_{i=18}^{33}V_{i}\leq 0.4\bullet\sum_{i=1}^{46}V_{i}\)(51)
\(0.09\leq V_{48}\leq 0.13\) (52)

Iterative Model Adoption and Optimization Solution Strategy

The identified rigorous mechanistic models for \(LD_{50}\), flash point, and odor type, the short-cut model for transparency, the surrogate model for sensorial rating as well as the heuristics in Eq. 46-52 are integrated to form the perfume formulation problem below.
\(\operatorname{}q_{s}\) (53)
s.t. Eq. 21-23 ANN-based surrogate model for \(q_{s}\)
Eq. 16-18 Design targets
Eq. 24-45 Mechanistic models
Eq. 46-52 Heuristics
Eq. 19-20 Design variables
This problem is implemented in GAMS 24.7 on a laptop with Intel 3.30 GHz CPU. The global solver BARON is used first and then the local solver SBB is employed if no optimal solutions are obtained from BARON.

GDP reformulation

Because of the complexity of the identified models and the number of intermediate variables, the problem is directly programmed using GDP. The disjunction is explicitly expressed as
\(\par \begin{bmatrix}Y_{i}\\ \par \begin{matrix}VL_{i}\leq V_{i}\leq VU_{i}\\ \par \begin{matrix}Eq.25\\ \par \begin{matrix}Eq.27-31\\ Eq.37-44\\ \end{matrix}\\ \end{matrix}\\ \end{matrix}\\ \end{bmatrix}\bigvee\par \begin{bmatrix}\neg Y_{i}\\ \par \begin{matrix}V_{i}=0\\ \par \begin{matrix}m_{i}=0\\ \par \begin{matrix}\text{FPLF}L_{i},FPP_{i},FP\gamma_{i},x_{i},FPPsat_{i}=0\\ n_{i,t},y_{i,t,z},\gamma_{i,t},x_{i,t},c_{i},\psi_{i}=0\\ \end{matrix}\\ \end{matrix}\\ \end{matrix}\\ \end{bmatrix}\) (54)
The GDP problem is further reformulated using the big-M approach with the solver JAMS and then solved by SBB. Different initial guesses are utilized. The second column of Table 5 lists the computational statistics. It contains 46 discrete variables, 9783 single variables, and 18230 equations. It takes 3459 seconds to obtain a local optimal solution.