Introduction

Most, if not all, manufactured products must meet quality specifications. These specifications can be roughly split into the following categories:

Legal : Imposed by governments in order to regulate a specific product/industry.

: Imposed by governments in order to regulate a specific product/industry. Technical : These requirements help us make sure that the product will work for the specific application it was designed.

: These requirements help us make sure that the product will work for the specific application it was designed. Commercial: When products are custom tailored for a specific customer, special quality demands can exist.

If the quality specification of a product are not met, the manufacturer incurs in a high economic penalty because the product cannot be sold, probably leading to a stock-out scenario. On the other hand, if the quality targets are overshot the manufacturer also incurs in an economic penalty (although not as high as in the previous scenario) because customers generally do not pay extra for that quality bonus. This is generally called quality giveaway.

In order to maximize profit, we need to find the process variables that balance the economic penalty of quality giveaway with the risk of manufacturing products below the minimum standard. This can sometimes be a difficult task because manufacturing processes can have many uncertainty sources.



Figure 1: Being able to quantify the uncertainty of decision variables allow us to optimize the profits of the manufacturing process.﻿

In the following example, we will use the Python package PyMC3 to build a probabilistic model of the blending process of gasoline and try to predict the amount of benzene in the feedstocks. The results of our model will yield useful information for the blending scheduler that will allow him/her keep quality giveaway to a minimum while also minimizing the risk of not meeting product specifications.

In general, the advantages of applying this type of algorithms include:

Automatic calibration of models to real world data (applicable not just to engineering or scientific applications).

Maximizing the utility of available information.

Understanding the variability of a process and identifying improvement opportunities.

Performing risk-based decision making.

Gasoline Blending

The gasoline blending process has the following steps:

A receipt is made by the blending scheduler. This receipt has the information of how much of every component needs to be pumped into the preparation tank. Each preparation tank has some residual content from the previous preparation that can take up to 20% of the total tank volume. There exists a maximum legal specification of 1% in volume of benzene for gasoline.

Each component is pumped separately from its tank into the preparation tank. Components have different properties that are characteristic of the manufacturing process that produced them. These properties can vary greatly from process to process but should stay approximately constant for a given process. Since component tanks are not regularly sampled and analyzed, the real benzene content is not known but because we know how that component was manufactured, we can make an educated guess. These are the properties we would like to estimate with our probabilistic model.

The tank is thoroughly mixed.

A sample is taken from the tank and analyzed in the lab following a very standardized procedure. The uncertainty regarding this analysis is very well understood and quantified: the repeatability is 0.1% in volume of benzene.

If the tank meets quality standards, some of it is sold and the remaining is blended into the following batch as described in step 1.

The Mathematical Model

For this particular project, we will follow sequential preparations of gasoline into an unique tank. We are going to assume the following hypotheses for building our mathematical model:

Benzene content of a mixture is linear with volume. This means that we can use the following equation to predict the benzene content of our blending tank:

where is the mixture benzene volume content, is the benzene volume content of the i-th component and is the volume added of the i-th component. The subscript 0 denotes the volume and properties of the initial volume of liquid of the tank.

There is no uncertainty associated to the pumping of the components to the blending tank: we assume that we know for certain how much of each component was pumped into the blending tank.

The lab result of the tank prepared at time t determines the amount of benzene in the residual contents of the tank prepared at time t+1 (this is how we model the sequential nature of the system).

Based on experience of the manufacturing process, resulting feed-stock quality stay pretty much constant within a three days time window. Because of this, we will use preparations t to t+2 to predict the outcome of preparation t+3. By studying the quality of our predictions we could determine the optimal time window for running the algorithm but this is out of the scope of the current project.

Errors of the lab analysis are normally distributed with mean = 0 and a standard deviation of 0.03 (this means that 99% of the times, results of the same sample should be at most 0.1% in volume apart).

We know how the manufacturing processes of the components used work and as a result we know how much benzene to expect in each tank. This will allow us to specify a prior probability distribution for our blending components. These variables are going to be described by Beta distributions (this distribution is specially useful for modelling bounded variables). Continuing with the notation so far:

where and are shape parameters that allow us to assign the right probabilities to the different possible values.



Figure 2: Shapes of the beta probability density function with different shape parameters.﻿

Based on the previously mentioned hypotheses, we can write the matrix form of the equation describing the behavior of the blending process as follows:

where is an by indexing the volumes of the m-th component in the preparation at the t-th time step( denotes the time period for which we consider that the properties of the blending components should stay the same), is a column vector indexing the volumetric benzene concentration for the m-th component, is a column vector of elements containing the volume of residual gasoline at the beginning of the blending process, is a column vector of elements with the volumetric benzene content of those residual volumes and is a column vector of elements containing the observed benzene content of the mixtures yielded by the laboratory analysis. The operator denotes the vector element-wise product.

Because of the residual content of the tanks, the final lab results for a given mixture becomes the starting point for the next one. We can express this relationship by defining the column vector , which contains all the laboratory analyses at time :

then,

One of the hypotheses previously enumerated was that the error of the lab analyses was normally distributed with and . This allows us to express every as

We can then write the final form of our blending equation that is going to be implemented in the software later on:

where:

Solving the Inference Problem

Imaging the following experiment:

You give a friend two fair six faced dice. You ask him to roll the dice and report the value of the sum of both dice without telling nor hinting their individual values. Your task is to make the best guess possible of what values each dice yielded based on the value of their sum.

Let’s imagine for a moment that your friend says he rolled a total of 8. What are the possible values for the dice? Any person would tell you 6+2, 5+3, 4+4, 3+5 and 2+6. What about the solutions 7+1 or -1 + 9? The latter set of results, although mathematically feasible, would be impossible because a six faced dice could never roll a -1 nor a 7. We could encode this information about the dice by using a uniform distribution that assigns a probability of 0 to all values that are impossible.



Figure 3: Probabilities of rolling a certain value for a fair six faced dice.

Going back to the set of feasible solutions, if the dice are perfectly fair, it is impossible for us to tell whether we have rolled a 2+6 or a 5+3. All the possible solutions have the same probability of happening.

Imagine that we switch one of the dice for another that has a bias for rolling a 3. The probability distributions for this new scenario would look like this:



Figure 4: Probabilities of rolling different values for the biased and fair dice.

Now, we know that it is much more likely for the orange dice to roll a 3, so if the total sum of both dice is 8, we could infer that the most likely solution is 5 + 3. Other solutions are still viable but probably less likely: rolling a 2 with the orange dice is less likely than rolling a 3 but still more likely than rolling a 5. All of this allows us to build a list of all the possible mathematical solutions to the problem and rank them based on the probability of happening.

For solving our blending model, we will use an algorithm from a family of algorithms called Markov Chain Monte Carlo. This will allow us to roll multiple dice for the variables that have uncertainty and carry out the simulations of the mixing process. The bias of each dice is going to be represented by our prior knowledge of the manufacturing process (the different beta distributions). After each roll, we will take not of the resulting benzene content of the mixture and compare it to the real lab result. If both of them agree within a certain margin (given by the repeatibility of the lab protocol), we will save the values of the input variables that produced such successful simulation. The more we save a given value, the higher ranking it will have and, as a result, the probability of it being the true magnitude of the variable will increase.

Implementation

The following code excerpt demonstrate how to implement a simpler version of the model described in the previous section. The code published in this article is not supposed to work nor to replicate the results shown in the figures and it is only for demonstrating a possible solution to the problem described above.

For this model, we are going to simulate a blending process consisting of four tanks containing gasoline from a catalytic cracker, light and heavy platformer naphtha and isomerate.

The first step is to define a way to obtain the appropriate beta distribution for specifying the input priors. We are going to crate a function that takes a minimum, maximum and average value together with the length of the expected range of values that our variable could take and returns the shape parameters of the beta distribution.

def get_beta_params(x_min, x_max, x_mu, x_range): var_loc = x_min var_scale = x_max-x_min var_mu = (x_mu-x_min)/var_scale var_var = ((x_range/var_scale)/3.0)**2 a = ((1-var_mu)/var_var - 1/var_mu)*var_mu**2 b = a * ((1/var_mu) - 1) return a, b, var_loc, var_scale

We can use the previous function as follows:

from scipy.stats import beta import numpy as np import matplotlib.pyplot as plt #Get the shape parameters a, b, loc, scale = get_beta_params(0.0, 10.0, 5.0, 3.0) x = np.arange(0, 1, 0.01) #Build the function's domain y = loc + scale * beta.pdf(x, a, b) #Get the probability density function #Plot plt.plot(x, y) plt.show()



Figure 5: Beta function produced by using the coefficients for the get_beta_params functions.﻿

The next step is to get the data for running the model.

#Dictionary containing the minimum, maximum, average and range values for the beta distributions of each priors tanks = { 'ccu': (0.0, 3.0, 1.0, 0.5), #Catalytic cracker gasoline 'cbo': (...), #Light platforming naphtha 'ptf': (...), #Heavy platforming naphtha 'iso':(None, None, None, None), #Isomerate. It has no benzene. } #Matrix indexing volumetric percentage used for each preparation per row #and each component per column. #Most recent preparations are going to be at the top of the rows. V = [] #Column vector with the starting volume percent of each preparation. #Most recent preparations are going to be at the top of the rows. Vb= [] #Column vector with the benzene of the initial tank's volume. #Most recent preparations are going to be at the top of the rows. xb= [] #Column vector with the benzene of the starting volumes of each preparation. #Most recent preparations are going to be at the top of the rows. lab_results = []

The next step is to use PyMC3 to build the probabilistic model and carry out inference to obtain the posterior distribution for our inputs.

import pymc3 as pm from Theano import tensor as T model = pm.Model() n_obs = len(b) with model: components_benzene = [] for t in tanks: if t == 'iso': x = 0.0 #Isomarate has no benzene. else: params = tanks[t] #Get shape parameters. a, b, loc, scale = get_beta_params(params[0], params[1], params[2], params[]) #Benzene content of the tank var = pm.Beta(t+'_var', alpha=a, beta=b, testval=0.5) x = pm.Deterministic(t, loc + scale * var) #Scaling to maximum and minimum values components_benzene.append([x]) x = T.stack(prop_list) # Column vector containing all the benzene variables # Blending simulation V_d = T.dot(V, x) # Linear blend b_d = Vb * xb # Bottom's contribution to blend y = pm.Deterministic('sim_data', V_d + b_d) # Blend all #Use lab data to infer posterior distributions likelihood = pm.Normal('benzene', mu=y, sd=0.05, observed=lab_results) start = pm.find_MAP() #Find a good estimate for the most probable outcome #Using the start value as a seed for MCMC trace = pm.sample(3000, start=start) #Infer

The main difference between this model and the one used for obtaining the results is that the latter was general enough to handle any number of components tanks with different types of products. Data was obtained by querying the laboratory and tank movements’ databases for records associated with every preparation. Once the data was collected, it was sequentially fed into the algorithm to obtain the results (making use of Theano’s shared variables allowed for significant speed ups of the code).

Results

The following figures show the results of the predictions using the naive approach of choosing the maximum likelihood values for each variables:



Figure 6: Prediction results made by the algorithm using the maximum likelihood approach. This naive approach yields predictions that are mostly within an error window of 0.15% in volume of benzene.﻿

It can be seen that the prediction (dotted blue line) follows the lab results (green line) with residuals oscillating around zero and amplitudes around 0.15%. Considering the laboratory protocol for determining the volumetric benzene content of gasoline has a repeatibility of 0.1%, the results obtained are extremely accurate.

Further investigation into the results shown in figure 4 showed that some of the volumes used for the estimation of the blending procedure where off by approximately 20% (specially in the amount of light platforming naphtha). This explains why there are some predictions that underestimate the benzene content by more than 0.2%. This could have been prevented by incorporating the uncertainty of the pumped volumes into the model or by using a different data source (tank radar data instead of accumulated flow meter readings) for getting the volume content of each component in the different preparations.

Next steps

If the economic penalty of having a stock-out and the quality giveaway is known, the results can be improved by implementing a loss function (also known as a cost function). The loss function would allow us to determine the optimal risk as a function of the profit given the probability distribution function that describes the uncertainty of our input variables.



Figure 7: Example of a cost function. These functions can be specified in a piece-wise manner in order to reflect the difference between the different scenarios.﻿

Another improvement that could be easily implemented would be updating the prior distribution of our input variables based on the results of the previous runs of the algorithm. This will allow the program to learn from previous preparations and improve the quality of predictions.

Finally, the algorithm showcased in this post can be easily generalized for predicting other properties, even when the mixing rules are non-linear.

Final Thoughts

This article shows a project that serves as a proof of concept for using probabilistic programming and Bayesian inference as tools for quantifying uncertainty in manufacturing processes. The results of applying this methodology to a real world problem show a promising outlook for this kind of algorithms and prove how relevant they can become in the near future.

If you liked the post or have any questions, please leave a comment!