# Performing an optimization in Python

Sun 13 June 2021This 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:

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.

# 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.