Add example of fitting binary model

Closes #71

doc/source/fitting/MultifluidBIPFitting.ipynb

@@ -0,0 +1,196 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "2dab4a3e",
+   "metadata": {},
+   "source": [
+    "# Multi-fluid Parameter Fitting\n",
+    "\n",
+    "Here is an example of fitting the $\\beta_T$ and $\\gamma_T$ values for the binary pair of propane+$n$-dodecane with the multi-fluid model. It uses differential evolution to do the global optimization, which is probably overkill in this case as the problem is 2D and other algorithms like Nelder-Mead or even approximate Hessian methods would probably be fine.\n",
+    "\n",
+    "In any case, it takes a few seconds to run (when the actual optimization is uncommented), demonstrating how one can fit model parameters with existing tooling from the scientific python stack."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3f8d8bce",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import json \n",
+    "import teqp, numpy as np, pandas, matplotlib.pyplot as plt\n",
+    "import scipy.interpolate, scipy.optimize\n",
+    "\n",
+    "import pandas\n",
+    "data = pandas.read_csv('VLE_data_propane_dodecane.csv')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b97637e7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def cost_function(parameters:np.ndarray, plot:bool=False):\n",
+    "\n",
+    "    # Fitting some parameters and fixing the others\n",
+    "    betaV, gammaV = 1.0, 1.0\n",
+    "    betaT, gammaT = parameters\n",
+    "\n",
+    "    # betaT, gammaT, betaV, gammaV = parameters\n",
+    "\n",
+    "    BIP = [{\n",
+    "        'function': '',\n",
+    "        'BibTeX': 'thiswork',\n",
+    "        'CAS1': '112-40-3',\n",
+    "        'CAS2': '74-98-6',\n",
+    "        'F': 0.0,\n",
+    "        'Name1': 'n-Dodecane',\n",
+    "        'Name2': 'n-Propane',\n",
+    "        'betaT': betaT,\n",
+    "        'betaV': betaV,\n",
+    "        'gammaT': gammaT,\n",
+    "        'gammaV': gammaV\n",
+    "    }]\n",
+    "    model = teqp.build_multifluid_model([\"n-Dodecane\", \"n-Propane\"], teqp.get_datapath(), \n",
+    "        BIPcollectionpath=json.dumps(BIP)\n",
+    "    )\n",
+    "    ancs = [model.build_ancillaries(ipure) for ipure in [0,1]]\n",
+    "\n",
+    "    cost = 0.0\n",
+    "    \n",
+    "    # The 0-based index of the fluid to start from. At this temperature, only one fluid \n",
+    "    # is subcritical, so it has to be that one, but in general you could start \n",
+    "    # from either one.\n",
+    "    ipure = 0 \n",
+    "\n",
+    "    for T in [419.15, 457.65]:\n",
+    "        # Subset the experimental data to match the isotherm \n",
+    "        # being fitted\n",
+    "        dfT = data[np.abs(data['T / K90'] - T) < 1e-3]\n",
+    "\n",
+    "        if plot:\n",
+    "            plt.plot(1-dfT['x[0] / mole frac.'], dfT['p / Pa']/1e6, 'X')\n",
+    "            plt.plot(1-dfT['y[0] / mole frac.'], dfT['p / Pa']/1e6, 'X')\n",
+    "\n",
+    "        try:\n",
+    "            # Get the molar concentrations of the pure fluid\n",
+    "            # at the starting point\n",
+    "            anc = ancs[ipure]\n",
+    "            rhoL0 = np.array([0, 0.0])\n",
+    "            rhoV0 = np.array([0, 0.0])\n",
+    "            rhoL0[ipure] = anc.rhoL(T)\n",
+    "            rhoV0[ipure] = anc.rhoV(T)\n",
+    "\n",
+    "            # Now we do the trace and convert retuned JSON\n",
+    "            # data into a DataFrame\n",
+    "            df = pandas.DataFrame(model.trace_VLE_isotherm_binary(T, rhoL0, rhoV0))\n",
+    "            \n",
+    "            if plot:\n",
+    "                plt.plot(df['xL_0 / mole frac.'], df['pL / Pa']/1e6)\n",
+    "                plt.plot(df['xV_0 / mole frac.'], df['pL / Pa']/1e6)\n",
+    "\n",
+    "            # Interpolate trace at experimental pressures along this \n",
+    "            # isotherm to get composition from the current model\n",
+    "            # The interpolators are set up to put in NaN for out\n",
+    "            # of range values\n",
+    "            x_interpolator = scipy.interpolate.interp1d(\n",
+    "                df['pL / Pa'], df['xL_0 / mole frac.'], \n",
+    "                fill_value=np.nan, bounds_error=False\n",
+    "            )\n",
+    "            y_interpolator = scipy.interpolate.interp1d(\n",
+    "                df['pL / Pa'], df['xV_0 / mole frac.'], \n",
+    "                fill_value=np.nan, bounds_error=False\n",
+    "            )\n",
+    "            # The interpolated values for the compositions \n",
+    "            # along the trace at experimental pressures\n",
+    "            x_model = x_interpolator(dfT['p / Pa'])\n",
+    "            y_model = y_interpolator(dfT['p / Pa'])\n",
+    "            if plot:\n",
+    "                plt.plot(x_model, dfT['p / Pa']/1e6, '.')\n",
+    "            \n",
+    "            # print(x_model, (1-dfT['x[0] (-)']))\n",
+    "            \n",
+    "            errTx = np.sum(np.abs(x_model-(1-dfT['x[0] / mole frac.'])))\n",
+    "            errTy = np.sum(np.abs(y_model-(1-dfT['y[0] / mole frac.'])))\n",
+    "            \n",
+    "            # If any point *cannot* be interpolated, throw out the model,\n",
+    "            # returning a large cost function value.\n",
+    "            #\n",
+    "            # Note: you might need to be more careful here,\n",
+    "            # if the points are close to the critical point, a good model might\n",
+    "            # (but not usually), undershoot the critical point of the \n",
+    "            # real mixture\n",
+    "            # \n",
+    "            # Also watch out for values of compositons in the data that are placeholders\n",
+    "            # with a value of nan, which will pollute the error calculation\n",
+    "            if not np.isfinite(errTx):\n",
+    "                return 1e6\n",
+    "            if not np.isfinite(errTy):\n",
+    "                return 1e6\n",
+    "            cost += errTx + errTy\n",
+    "\n",
+    "        except BaseException as BE:\n",
+    "            print(BE)\n",
+    "            pass \n",
+    "    if plot:\n",
+    "        plt.title(f'dodecane(1) + propane(2)')\n",
+    "        plt.xlabel('$x_1$ / mole frac.'); plt.ylabel('$p$ / MPa')\n",
+    "        plt.savefig('n-Dodecane+propane.pdf')\n",
+    "        plt.show()\n",
+    "\n",
+    "    return cost"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4ee30326",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The final parameter values, will be overwritten if \n",
+    "# optimization call is uncommented\n",
+    "x = [1.01778992, 1.17318854]\n",
+    "\n",
+    "# Here is the code used to do the optimization, uncomment to run it\n",
+    "# Note: it is commented out because it takes too long to run on doc builder\n",
+    "#\n",
+    "# res = scipy.optimize.differential_evolution(\n",
+    "#     cost_function, \n",
+    "#     bounds=((0.9, 1.5), (0.75, 1.5)), \n",
+    "#     disp=True, \n",
+    "#     polish=False\n",
+    "# )\n",
+    "# print(res)\n",
+    "# x = res.x\n",
+    "\n",
+    "cost_function(x, plot=True)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

doc/source/fitting/VLE_data_propane_dodecane.csv

@@ -0,0 +1,27 @@
+Ncomp,T / K90,fluid[0],fluid[1],fluid[2],p / Pa,x[0] / mole frac.,x[1] / mole frac.,x[2] / mole frac.,y[0] / mole frac.,y[1] / mole frac.,y[2] / mole frac.,p / MPa
+2,419.15,PROPANE,C12,,453000.000,0.0755,0.9245,,0.94865,0.05135,,0.453
+2,419.15,PROPANE,C12,,622000.000,0.1048,0.8952,,0.95308,0.04692,,0.622
+2,419.15,PROPANE,C12,,1011000.000,0.1642,0.8358,,0.96513,0.03487,,1.011
+2,419.15,PROPANE,C12,,1694000.000,0.2766,0.7234,,0.97512,0.02488,,1.694
+2,419.15,PROPANE,C12,,2185000.000,0.3403,0.6597,,0.97512,0.02488,,2.185
+2,419.15,PROPANE,C12,,2908000.000,0.4651,0.5349,,0.97518,0.02482,,2.908
+2,419.15,PROPANE,C12,,3397000.000,0.5334,0.4666,,0.97598,0.02402,,3.397
+2,419.15,PROPANE,C12,,3907000.000,0.5938,0.4062,,0.97276,0.02724,,3.907
+2,419.15,PROPANE,C12,,4820000.000,0.7033,0.2967,,0.97255,0.02745,,4.820
+2,457.65,PROPANE,C12,,406000.000,0.0451,0.9549,,0.8682,0.1318,,0.406
+2,457.65,PROPANE,C12,,663000.000,0.0872,0.9128,,0.90189,0.09811,,0.663
+2,457.65,PROPANE,C12,,830000.000,0.0999,0.9001,,0.91889,0.08111,,0.830
+2,457.65,PROPANE,C12,,1367000.000,0.1698,0.8302,,0.94423,0.05577,,1.367
+2,457.65,PROPANE,C12,,1825000.000,0.2306,0.7694,,0.95299,0.04701,,1.825
+2,457.65,PROPANE,C12,,2309000.000,0.2901,0.7099,,0.95408,0.04592,,2.309
+2,457.65,PROPANE,C12,,2815000.000,0.3293,0.6707,,0.96007,0.03993,,2.815
+2,457.65,PROPANE,C12,,3813000.000,0.4361,0.5639,,0.96134,0.03866,,3.813
+2,457.65,PROPANE,C12,,4792000.000,0.5363,0.4637,,0.95919,0.04081,,4.792
+2,457.65,PROPANE,C12,,5799000.000,0.6278,0.3722,,0.95279,0.04721,,5.799
+2,457.65,PROPANE,C12,,6640000.000,0.7008,0.2992,,0.94045,0.05955,,6.640
+,,,,,,,,,,,,
+,,,,,,,,,,,,
+,,,,,,,,,,,,
+,,,,,,,,,,,,
+,,,,,,,,,,,,
+,,,,,,,,,,,,

doc/source/fitting/index.rst

@@ -0,0 +1,9 @@
+Fitting
+=======
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+   MultifluidBIPFitting
+

doc/source/index.rst

@@ -11,6 +11,7 @@ Welcome to teqp's documentation!
    derivs/index
    algorithms/index
    examples/index
+   fitting/index
    api/modules
 
 Indices and tables