Add example of parameter fitting with MF to teqp docs

[ianhbell]

For instance two isotherms for a binary mixture, something like propane + octane or something like that.

[ianhbell]

import json 
import teqp, numpy as np, pandas, matplotlib.pyplot as plt
import scipy.interpolate, scipy.optimize

import pandas
data = pandas.read_excel('VLE_data_propane_dodecane.xlsx')

def cost_function(parameters:np.ndarray, plot:bool=False):

    # Fitting some parameters and fixing the others
    betaV, gammaV = 1.0, 1.0
    betaT, gammaT = parameters

    # betaT, gammaT, betaV, gammaV = parameters

    BIP = [{
        'function': '',
        'BibTeX': 'thiswork',
        'CAS1': '112-40-3',
        'CAS2': '74-98-6',
        'F': 0.0,
        'Name1': 'n-Dodecane',
        'Name2': 'n-Propane',
        'betaT': betaT,
        'betaV': betaV,
        'gammaT': gammaT,
        'gammaV': gammaV
    }]
    model = teqp.build_multifluid_model(["n-Dodecane", "n-Propane"], teqp.get_datapath(), 
        BIPcollectionpath=json.dumps(BIP)
    )
    ancs = [model.build_ancillaries(ipure) for ipure in [0,1]]

    cost = 0.0
    
    # The 0-based index of the fluid to start from. At this temperature, only one fluid 
    # is subcritical, so it has to be that one, but in general you could start 
    # from either one.
    ipure = 0 

    for T in [419.15, 457.65]:
        # Subset the experimental data to match the isotherm 
        # being fitted
        dfT = data[np.abs(data['T (K90)'] - T) < 1e-3]

        if plot:
            plt.plot(1-dfT['x[0] (-)'], dfT['p (Pa)']/1e6, 'X')
            plt.plot(1-dfT['y[0] (-)'], dfT['p (Pa)']/1e6, 'X')

        try:
            # Get the molar concentrations of the pure fluid
            # at the starting point
            anc = ancs[ipure]
            rhoL0 = np.array([0, 0.0])
            rhoV0 = np.array([0, 0.0])
            rhoL0[ipure] = anc.rhoL(T)
            rhoV0[ipure] = anc.rhoV(T)

            # Now we do the trace and convert retuned JSON
            # data into a DataFrame
            df = pandas.DataFrame(model.trace_VLE_isotherm_binary(T, rhoL0, rhoV0))
            
            if plot:
                plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6)
                plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6)

            # Interpolate trace at experimental pressures along this 
            # isotherm to get composition from the current model
            # The interpolators are set up to put in NaN for out
            # of range values
            x_interpolator = scipy.interpolate.interp1d(
                df['pL / Pa'], df['xL_0 / mole frac.'], 
                fill_value=np.nan, bounds_error=False
            )
            y_interpolator = scipy.interpolate.interp1d(
                df['pL / Pa'], df['xV_0 / mole frac.'], 
                fill_value=np.nan, bounds_error=False
            )
            # The interpolated values for the compositions 
            # along the trace at experimental pressures
            x_model = x_interpolator(dfT['p (Pa)'])
            y_model = y_interpolator(dfT['p (Pa)'])
            if plot:
                plt.plot(x_model, dfT['p (Pa)']/1e6, '.')
            
            # print(x_model, (1-dfT['x[0] (-)']))
            
            errTx = np.sum(np.abs(x_model-(1-dfT['x[0] (-)'])))
            errTy = np.sum(np.abs(y_model-(1-dfT['y[0] (-)'])))
            
            # If any point *cannot* be interpolated, throw out the model,
            # returning a large cost function value.
            #
            # Note: you might need to be more careful here,
            # if the points are close to the critical point, a good model might
            # (but not usually), undershoot the critical point of the 
            # real mixture
            # 
            # Also watch out for values of compositons in the data that are placeholders
            # with a value of nan, which will pollute the error calculation
            if not np.isfinite(errTx):
                return 1e6
            if not np.isfinite(errTy):
                return 1e6
            cost += errTx + errTy

        except BaseException as BE:
            print(BE)
            pass 
    if plot:
        plt.title(f'dodecane(1) + propane(2)')
        plt.xlabel('$x_1$ / mole frac.'); plt.ylabel('p / MPa')
        plt.savefig('n-Dodecane+propane.pdf')
        plt.show()
    return cost

if __name__ == '__main__':
    res = scipy.optimize.differential_evolution(cost_function, bounds=((0.9, 1.5), (0.75, 1.5)), disp=True, polish=False)
    print(res)
    cost_function(res.x, plot=True)