The optimal solution of a mixture problem with modeFRONTIER

Table 1: The table collects some examples of engineering fields where mixture design problems could arise.

A simple mixture problem
As explained before, in a simple mixture design problem the final result does not depend on the quantity of the ingredients but rather on their proportions. This means that the problem can be formulated in terms of relative concentrations ci of the ingredients, which are subjected to the following constraint:

where n is the number of the ingredients in the mixture and qi is the quantity of the ith ingredient. It is clear that the violation of this constraint leads to a physical non-sense, hence it is mandatory that the solution of a mixture problem
satisfies this constraint; to this aim, it is recommended to rewrite the above equation in the following form, in order to facilitate the problem solution.

In other words, one ingredient concentration (e.g. the last) is not a free parameter but it can be computed once the others are known. It always has to belong to the interval [0,1] to maintain a physical meaning.
The mixture design problem can be described schematically as drawn in Figure 1, where, on the left, the unknown concentrations and, on the right, the results of the mixture are displayed; the central box indicates that a process (chemical, physical or whatever makes sense in the analyzed context) transforms the inputs into outputs.
As mentioned above, the solution of such problems is usually more challenging in the presence of many objectives and constraints involving the mixture results. They can be in contrast to one another and the constraints can seriously restrict the possibilities to find out a proper solution.
A mixture design, as described before, can be regarded as an optimization problem where the ingredients' concentrations are the inputs, the mixture results are the outputs and all the requisites can be transformed into objectives (maximizations or minimizations) and constraints.
In this work the following simple mixture problem has been tackled: five ingredients (let us say A - red, B - green, C- blue, D – deep blue and E - black) have to be mixed together to obtain, within a given time, a solution characterized by the highest concentration possible for E. The ingredients involve different costs (this is quite typical) and therefore, we are also looking for the most economic recipe.

Figure 1: A schematic representation of a simple mixture problem. On the left the unknown ingredients’ concentrations are displayed. On the right the products of a mixture process are represented together with the requisites (objectives and constraints) that the mixture should satisfy.

Moreover, the components B and E can be extracted from the final mixture; E will be the product of the mixture process we are looking for, while B can be considered a sub-product and it can be reused in a new production cycle. The other products (if any) cannot be reused and they have to be considered as a discard.

In Table 2 the possible reactions between the ingredients and the costs for a unit quantity of the ingredients are listed.

Table 2: The table collects the possible reactions that can take place when two particles reach the same position and the ingredient costs. For example, if a particle B and a particle D reach the same position, they will be transformed into two particles of E. Ingredient E is considered as the main result of the mixture process while B can be extracted from the mixture and reused for a new production cycle.

Another issue has to be considered: as the mixing process proceeds, the probability that the reactions described in Table 2 take place decreases up to zero with a quadratic law, as shown in Figure 2. This means that the process tends to a
stable point which corresponds to assume an invariable status after a certain time (which we will call "reaction time" from now on).

Figure 2: The probability that the reactions between the ingredients described in Table 2 follows the law shown in the picture; a first unitary constant plateau is followed by a quadratic decreasing piece, up to the reaction time. It is clear that other decay laws could be easily adopted.

The mixture process is simulated in this way: an array of 250 x 250 cells has been considered as a sort of terrain where the mixture process can take place. The ingredients are modeled by means of a given number of particles (computed according to the concentration) which are randomly positioned in the cells array at the beginning of the mixture process (let us say at time zero). The initialization does not allow that more than one particle falls into a single cell. The cells array is not completely filled by the particles, only 25% of it will be occupied during the simulation; this choice is obviously arbitrary and it could be easily changed if necessary.
Once the initialization has been concluded, time starts, it increases its value by a unit at a time, and all the particles in the array are randomly moved in one of the eight possible surrounding cells to obtain a new configuration.

It obviously could happen that two particles fall in the same cell: In this case, a reaction could take place. A random number in the interval [0,1] is generated and if it is less than the actual reaction probability (see Figure 2), the reaction starts according to the rules summarized in Table 2 involving two particles at a time. When more than two particles reach the same cell the reactions are performed sequentially starting from the first couple of particles up to the last one.

Figure 3: The curve fitting tool has been used to show how the final concentrations of ingredients are distributed. In the picture the concentration of the ingredient E is reported, being the following initial concentrations: A = 0.1, B = 0.3, C = 0.4, D = 0.2 and E = 0 for the first case (left) and A = 0.0, B = 0.2, C = 0.5, D = 0.2 and E = 0.1 for the second case (right). It comes out that in the first case the Gaussian distribution has the highest Kolmogorov-Smirnov test score (83%) while in the second case the logistic distribution seems to be the best one (57%).

The mixture process variability
As mentioned above the initialization of the array is done using a random criterion: moreover, the particle motion is not given by a deterministic law, but it looks like a random walk in the array. To conclude, also the fact that a reaction can or cannot take place is governed by chance. This means that a mixture process cannot be exactly reproduced; if two simulations of the same mixture process (with the same initial concentrations) are performed, probably two different final results will be obtained.

When the process time has reached a given value (in our case 200) the mixture does not change anymore (no reaction can take place) and therefore the ingredient concentrations are computed simply by considering the number of ingredient particles in the array divided by the total number of particles. As the reader has probably understood by now, the reaction process has been simulated in such a way that the number of particles does not change during the mixture process, preserving in this way the total mass. The mixture process described above can be simulated using, for example, a Matlab script (other choices are obviously possible) or better a compiled routine (Fortran, C, …) to improve the performance. It is straightforward to note that the computational cost grows quadratically with the array dimensions, which have to be sufficiently large to enable a reliable simulation of the process, and this could lead to prohibitive solution times.

Figure 4: An example of the cells array during a mixture process at time t = 50 (particle colors are the same as used in Table 2). The following initial concentrations have been used: A = 0.1, B = 0.3, C = 0.4, D = 0.2, E = 0.

Figure 5: The ingredient concentrations plotted versus the process time. The following initial concentrations have been used: A = 0.1, B = 0.3, C = 0.4, D = 0.2 and E = 0.

These differences are clearly due to the sources of randomness which are present in the model, as mentioned above: however, their effects reduce as the simulation time increases and the mixture process tends in mean to a given final configuration; this is due to the fact that the reaction probability (see Figure 2) decreases up to zero with the process time, leading in this way to an inert mixture.

This makes everything more complicated; actually, an optimal solution could be characterized by important variations in the final concentrations which could be unacceptable.

At this point, it becomes mandatory to identify the probability density function (or the cumulative one) that characterizes the final concentrations, in order to correctly judge the goodness of the recipe. In this way, it is actually possible to estimate the final concentration of the ingredients with a probabilistic approach which appears to be, in this case, more reliable than a simple deterministic one. To do this, it is necessary to simulate more than once the mixture process and find out a statistical distribution that fits at best the results. This is a distribution fitting problem which is not exactly easy to solve in its more general formulation.

Table 3: The table collects the equations of the Gaussian and the logistic probability density functions (PDF) and the values of the location and scale parameter as estimated with the method of moments. The empirical mean and standard deviation are μ and σ respectively. The last column collects the equations that allow to compute the variable corresponding to a given value of cumulative distribution function P(x < ).

Once the parameters of both distributions have been computed, a Pearson chi-squared test can be performed to find out which distribution best fits the data.

Finally, the simulation of a given recipe has to include these steps; the mixture process simulation has to be repeated for a certain number of times and a dataset containing the final concentration for each ingredient has to be built. Then, the method of moments can be used to identify the probability density function parameters and a Pearson chi-squared test has to be performed to find out the best theoretical distribution, among the ones that have been considered. To conclude, the value that has a very low probability to overestimate the final ingredient concentration has to be computed.

All the steps described above have to be performed by the simulation software which will pass all the information pertaining to the process to modeFRONTIER.

In order to show that different recipes could lead, in general, to different final concentration distributions we have decided to consider these two different initial conditions; in the first case A = 0.1, B = 0.3, C = 0.4, D = 0.2 and E = 0 are the initial concentrations, while in the second we have A = 0.0, B = 0.2, C = 0.5, D = 0.2 and E = 0.1.
The mixture processes have been simulated by means of 100 runs for each recipe and the final concentrations have been stored on a text file; then, these data have been loaded in the Designs Space of modeFRONTIER and some statistics have been done.
As the Figure 3 shows, the histograms of the ingredient E are well fitted by the Gaussian, in the first case, and by the logistic distribution, in the second case; these distributions actually obtain the highest Kolmogorov-Smirnov scores (83% and 57% respectively).

The estimation of the mean and standard deviation of the final concentrations
As explained in the previous paragraph, the method of moments is used to determine the location and scale parameters of the Gaussian and the logistic distributions. Therefore, it becomes mandatory to perform a reliable estimation of the mean μ and the standard deviation σ of the final ingredient concentrations. The simplest approach is to use the definition of the estimated values of these quantities:

where ci is the concentration of a given ingredient while n is the number of runs which have been performed. The main drawback of this approach is that a large amount of runs n is usually needed to obtain a sufficiently good estimation, especially for the standard deviation, with the consequent deterioration of the computational time. For this reason the choice of an appropriate number of runs n is absolutely crucial.

In order to choose a proper number of runs the following approach has been used: it has been assumed that the estimations corresponding to 1000 runs give the right values of the mean and the standard deviation. Then, the difference between the mean and three times the standard deviation has been computed for all the estimations, supposing, for simplicity, that the output is always normally distributed (see Table 3); the obtained value represents the worst concentration for ingredient E we would get from the mixture process if the estimations of the mean and the standard deviation were correct. Actually, if we compute the absolute difference between such quantity and the analogous one evaluated with 1000 runs (let us call this quantity: the Error), we obtain an estimation of the absolute error we commit when computing the lowest concentration for the ingredient E we can expect from the mixture process.
Looking at Figure 6, where this Error is plotted against the number of runs, it is clear that after 100 runs the estimated error is always less than 1*10-3, which is the sensibility we have in this problem. For this reason and in view of the introduced approximations, 200 runs can be considered as a good choice for a sufficiently accurate estimation of the mean and the standard deviation for the ingredient E.

Figure 6: The Error, computed as the absolute difference between the mean and three times the standard deviation and the analogous quantity computed at 1000 runs. After 100 runs this quantity is always less than 1*10-3, which represents the sensibility on the mixture concentration measurements (A = 0.1, B = 0.3, C = 0.4, D = 0.2 and E = 0 have been set as initial concentrations).

Analogous results can be obtained considering different initial ingredient concentrations.

The modeFRONTIER solution
In order to efficiently solve the mixture design problem a workflow has to be defined following a similar structure of Figure 1.
We decided to use a step of 1*10-3 to discretize the initial concentrations. It is interesting to note that this step represents the smallest possible variation of a concentration which corresponds to a difference of 15 particles (on a total of 15625, with an array filling of 25%) in the cells array. Obviously, it does not make sense to use steps which correspond to a particle variation less that one and, however, it is reasonable to have at least a certain number of particles in the array, in order to have a model too sensitive to randomness.

As explained before, an ingredient concentration can be computed once the other concentrations are known, according to equation [2]. The resulting concentration has to fall in the interval [0,1] and therefore a constraint has to be added.
A calculator can be used to compute the cost of the mixture process taking into account the data in Table 2 and the fact that ingredients B and E can be reused for a future production cycle; to this aim a simple JavaScript can be used.

Two important constraints are given: the first one involves the production cost, which has to be less than 1. The second one states that the gain of ingredient E has to be greater than 0.2: we actually are not interested in a recipe, even a cheap one, which does not produce an interesting gain in ingredient E.

Figure 7: The modeFRONTIER workflow ready to run.

A random DOE of 20 initial mixtures is adopted to feed a MOGA-II optimization with 20 generations: however, one initial design has been modified to include the following initial condition: A = 0.0, B = 0.5, C = 0.0, D = 0.5 and E = 0.0 (design ID20), because, looking at Table 2, where the possible reactions between ingredients are reported, it turns out that a high initial concentration of B and D should produce a high final concentration for the ingredient E. A second design A = 0.0, B = 0.4, C = 0.2, D = 0.4 and E = 0.0 (design ID19) is introduced in the DOE table following the same consideration.

The default values in the Scheduler node for the algorithm set up are used. The unfeasible designs (the ones which do not satisfy equation [1]) are not evaluated, because they do not have any physical sense.

After the optimization process it is possible to find out the Pareto front solutions. In Figure 8 the production cost is plotted versus the minimum gain of ingredient E. The DOE designs are the blue points while the red points are the Pareto designs after 30 generations. It is interesting to note that the two initial designs inserted in the DOE table in the attempt to furnish good solutions are both unfeasible; they have a high final concentration of E but they are too expensive.

Moreover, there are some Pareto designs that have a negative cost and a gain in E greater than 0.2; surprisingly, this means that, using some recipes, one can organize a process able to produce a certain quantity of the ingredient E, according to the main objective, and have a production of B sufficient to cover all the other ingredients’ cost and more.

The optimal solutions are characterized by low or null values of initial concentration of ingredient E and relatively high of A. Initially, there was the suspect that ingredient B and D should have a high initial concentration; the optimization process gives evidence of the contrary, actually the optimal solutions usually have low initial concentrations of such ingredients.

The robustness of the optimal solution
Once the Pareto front has been found, thanks to the optimization process, it is necessary to choose just one configuration. Obviously, many different criteria can be used to choose the best-for-us configuration. Looking at Figure 8 it immediately becomes clear, that the cost of the process drastically increases when a gain in E greater than 0.29 is reached. For this reason we decided to adopt an optimal solution pertaining to the knee of the Pareto front which provides a reasonable high gain level without being too expensive. The choice falls to design number 462, which has the following initial concentrations: A = 0.769, B = 0.000, C = 0.209, D = 0.220 and E = 0.000. It can be further noted that this recipe has the interesting property to have only three initial ingredients with non-zero value concentrations; this is extremely interesting because it allows to enormously simplify the production process. The cost of such recipe (cost = 0.198) is definitely lower than the limit of 1 imposed to the optimization process and the production of ingredient E (gain_fE = 0.290) is one of the highest found.

Figure 8: The production cost plotted versus the minimum gain of ingredient E. The DOE designs are the blue points while the red points are the Pareto designs after 30 generations. All the designs (feasible and unfeasible) are reported here. The user defined DOE designs 18 and 19 are highlighted in green and they fall in the upper-right part of the scatter.

Another interesting issue that could come into play when designing a mixture is the robustness of the recipe. The production process could actually start with a slightly different recipe with respect to the optimal one; the reasons of these variations can arise, for example, from tolerances in machining, in measurements, in the ingredient physical properties and more.

The ingredient final concentrations could be affected in an unacceptable way by these uncertainties: the aim of this paragraph is mainly to show how it is possible to check if the optimal solution is or is not sensitive to the noise factors.
The first step is certainly to characterize the input variables, which represent the ingredients' initial concentrations, with probabilistic density functions that summarize the effects of all the noises.
Then, it is necessary to run a certain number of simulations of the optimal recipe taking into account the effects of noise, as imposed by the stochastic approach, and try to find out the probability density functions that best fit the outputs. Finally, the noise effects can be analyzed and the robustness of the recipe measured and eventually compared with restrictions or requirements that the process has to guarantee.

The problem is very similar to that exposed and solved before, with the important difference that now the recipe is known and that the initial concentrations are affected by small variations whose probability to appear is driven by some given PDFs.

The solution of such problem could be led as before; a certain number of simulations could be run and the final ingredient concentrations stored. Then, another fitting problem should be solved and the robustness of the solution estimated. As shown before, the most important drawback of this approach is that a large amount of repetitions is needed to have a good estimation of the mean and the standard deviation.
To partially mitigate this aspect one can use the polynomial chaos method (see [1]), which is known as a very accurate and relatively cheap approach.
We have supposed that the ingredients’ concentrations are all characterized by a Gaussian distribution with a standard deviation of 0.6*10-2; In modeFRONTIER it is easy to set up a robustness analysis; starting from the existing project the user has to change the input variable definitions to provide their correct stochastic definitions, copy in the DOE table the optimal configuration of which robustness has to be measured and finally set up the MORDO panel in the Scheduler node: this last step is probably the most important one, because the number of designs to be generated and evaluated around the nominal configuration and the polynomial chaos order have to be chosen. In our case we decided to generate 50 designs with a Latin-Hypercube technique and to use a chaos expansion of order 2.

Figure 9: The distribution fitting tool has been used to find out the probability density functions that better fit the mixture cost and gain of ingredient E. It can be seen that the Gaussian, the logistic and the Weibull distributions usually are the ones that have the highest Kolmogorov-Smirnov rates. Also gamma and beta distributions sometimes seem to fit the data in a proper way. Once plotted however, these functions are all very close to one another.

The distribution fitting tool can be used once again to find out the theoretical distribution that best fits the data pertaining to the final concentrations and to understand if the variations around the expected values are acceptable or not.

Table 4 collects some estimations of the mean and the standard deviation of the mixture cost and gain in ingredient E, which can be accessed using different modeFRONTIER tools. It immediately appears that the estimations of the means always provide the same values; unfortunately the same does not happen for the standard deviation. Very different numerical values can be obtained using different approaches. The polynomial chaos approach can be considered the most accurate one and it is strongly recommended to use this approach, when possible.

Table 4: The estimated means and standard deviations of the gain in E and the mixture cost have been estimated using three different approaches: in (1) equations [3] and [4] have been computed, in (2) the curve fitting tool has been used (the maximization of likelihood is used to determine the pdfs parameters: the normal distribution has been supposed to characterize both the gain and cost), while in (3) the polynomial chaos approach has been adopted. It can be seen that important differences in the estimated values of the standard deviation are present.

It can be further noted that the estimated values of the standard deviations computed by the polynomial chaos approach are greater than the analogous ones computed with the other techniques: these differences, which are not negligible in this case, could lead to an uncorrect estimation of the robustness of the recipe.

Actually, if we suppose that both the gain and the cost are normally distributed, it is possible to compute the worst scenario in terms of these two system outputs as summarized in Table 5. It can be easily seen that the worst scenario strongly depends on the technique (or, in other words, on the accuracy) adopted for the estimation of the mean and the standard deviation of the system outputs.

Table 5: The estimated worst scenarios for the Gain E and the mixture cost using different estimated values of the means and standard deviations. The two outputs have been both supposed to be normally distributed. It can be seen that the worst scenarios strongly depend on the technique (accuracy) used for the estimation of the mean and the standard deviation of the output distributions.

Conclusions
In this work a mixture design problem has been solved with modeFRONTIER; several aspects that can arise in a typical situation have been considered and some considerations on the process variability and solution robustness have been verified.

We have shown how complicated mixture problems can be tackled involving many ingredients, satisfying both physical and economical objectives and constraints. It is also possible to estimate the robustness of an optimal configuration and to understand in this way if the process noises can affect the final results in an undesired way.

References

• [1] Lovison Alberto (2008), Uncertainty Quantification in modeFRONTIER: Monte Carlo, Latin Hypercube Sampling and Polynomial Chaos, modeFRONTIER Technical Report
  
• [2] Stefan H. Steiner, Michael Hamada, Bethany J. Giddings White, Vadim Kutsyy, Sofia Mosesova, Geoffrey Salloum, (2007), A Bubble Mixture Experiment Project for Use in an Advanced Design of Experiments Class, Journal of Statistics Education, Volume 15, Number 1

Contacts
For more information on this article and topic, please contact the author: Massimiliano Margonari - Enginsoft S.p.A.
info@enginsoft.it