Performing an optimization in PythonSun 13 June 2021
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.
The Antoine equation is a semi-empirical correlation that describes the relation between vapor pressure and temperature for pure substances. It takes the form:
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 following plot presents the vapor pressure data for propane that was determined experimentally.
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
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)+B/(C+temp)-A eq2 = (A-np.log10(pvap))*(C+temp)-B eq3 = B/(A-np.log10(pvap))-temp-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
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-ABC/(ABC+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, ABC, ABC))
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-ABC/(ABC+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.