Performing an optimization in Python

This article covers a simple example of how to use Python to perform an optimization problem. We will fit the Antoine equation to experimental data in order to determine the Antoine constants of pure propane.

Antoine equation

The Antoine equation is a semi-empirical correlation that describes the relation between vapor pressure and temperature for pure substances. It takes the form:

$$\log(P^v) = A - \frac{B}{C + T}$$

where \(P^v\) is the vapor pressure in mmHg and \(T\) is the temperature in °C. \(A\), \(B\), and \(C\) are the Antoine constants and are specific for each substance.

The data

The following plot presents the vapor pressure data for propane that was determined experimentally.

Data

Fitting the Antoine equation

We will require a few Python libraries. numpy is the main scientific computing library for Python and scipy has the necessary optimization routines. We will also use matplotlib to create the plots.

# Imports
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

First, we place the data from the table above into a file named propane-vapor-pressure.txt, load it using numpy, and create temperature and pressure arrays.

# Load data
T_Pvap = np.loadtxt('propane-vapor-pressure.txt')
temp = T_Pvap[:,0]  # ºF
pvap = T_Pvap[:,1]  # psia

We will need to convert the data to the relevant units. Temperature from °F to °C to and vapor pressure from psia to mmHg.

# Unit conversion
temp = (temp-32)/1.8  # ºC
pvap = pvap/14.69595*760  # mmHg

In most numerical optimization routines we will have to provide initial estimates of the values we are trying to optimize.

In this example, we can find very reasonable estimates of the Antoine parameters A, B and C by solving the Antoine equation for any 3 known points. Here we will pick the 3 initial points and solve this system of 3 equations to obtain initial estimates. The obtained parameters will be true for those points but will have increasing deviations as we get further from these initial conditions.

To solve the resulting system of nonlinear equations we can use the fsolve solver from scipy. We will save our initial estimate in an array named ABC0.

Note that we could attempt to provide a random initial estimates, although this results in higher chances of failing the optimization. Good initial estimates can be very important, in some cases, to allow the optimization algorithm to converge to a relevant optimum.

# Define the system of 3 equations and solve it using scipy.optimize.fsolve
def equations(p):
    A, B, C = p
    eq1 = np.log10(pvap[0])+B/(C+temp[0])-A
    eq2 = (A-np.log10(pvap[1]))*(C+temp[1])-B
    eq3 = B/(A-np.log10(pvap[2]))-temp[2]-C
    return eq1, eq2, eq3

A, B, C = optimize.fsolve(equations, (5, 600, 200))

ABC0 = np.array([A, B, C])

Now that we have a good initial estimate we can run an optimization routine to find the set of parameters A, B and C that best fit our data. We use scipy fmin routine. This function will minimize the value of a given objective function.

For a given set of parameters, our objective function will calculate the vapor pressures for each temperature and compare those results to the experimental values we stored in the pvap array. We use a least square objective function here. fmin will then try to find the set of A, B and C that minimize this objective function.

# Define the objective function
def fobj(ABC):
    pvap_calc = 10**(ABC[0]-ABC[1]/(ABC[2]+temp))
    f = np.sum((pvap_calc-pvap)**2)
    return f

To run the optimization we simple call the fmin function and provide our objective function and the initial estimates. We will then print the results.

# Run optimization to determine parameters A, B and C
ABC = optimize.fmin(fobj, ABC0)

print('\nOptimized Antoine Constants:\nA = %.3f\nB = %.3f\nC = %.3f' % (ABC[0], ABC[1], ABC[2]))

If the optimization finished successfully you will see something like:

Optimization terminated successfully.
        Current function value: 1500.260802
        Iterations: 215
        Function evaluations: 391

Optimized Antoine Constants:
A = 7.002
B = 899.442
C = 260.441

Finally, we can calculate the values of vapor pressure that result from our optimized parameters and plot experimental and calculated results.

pvap_final = 10**(ABC[0]-ABC[1]/(ABC[2]+temp))

# Display results
fig, ax = plt.subplots()
ax.plot(temp, pvap, 'bo', label='Experimental')
ax.plot(temp, pvap_final, 'r-', label='Calculated')
plt.xlabel('Temperature (ºC)')
plt.ylabel('Vapor Pressure (mmHg)')
legend = ax.legend(loc='lower right')
plt.show()

Here is the plot we obtain. We achieve a very good fit.

Fit

Source code

Source code on Github

References

Wikipedia - Antoine equation