add tuning folder

degradation_model/tuning_code/readme.md

@@ -0,0 +1,43 @@
+# Code to tune degradation model
+
+The process to tune the degradation model is documented in this [paper](https://iopscience.iop.org/article/10.1149/1945-7111/ad1294) and is summarized in this figure:
+
+![degradation tuning summary](./tuning_summary.png)
+
+We adopt a sequential tuning process and it includes 4 steps:
+
+## Step 1
+1. Using calendar aging data for multiple temperatures to tune the SEI model. In this case, the contribution of mechanical damage and Li plating to aging would be negligible.
+2. The SEI model is tuned using the data from the calendar aged cells in both cold (-5°C) and hot temperature (45°C) conditions. The reason for doing this is to capture the temperature dependency of SEI growth, being faster at high temperatures and slower at colder temperatures.
+3. The vector of degradation mechanism parameters P to be tuned in step 1 are:
+$$P_\mathrm{cal} = \left[k_{0,\mathrm{SEI}},D_\mathrm{SEI},E_{a,\mathrm{SEI}}\right]^T $$
+4. The vector of data to be fitted against model outputs are:
+$$Y_\R\mathrm{cal} = \left[C, n_\mathrm{Li},y_0,x_{100}\right]^T$$
+
+To perform step 1, please run this [notebook](./step_1_calendar.ipynb)
+
+### Input Data Required:
+eSOH parameters $[x_0,x_{100},y_0,y_{100},C_n,C_p,C,n_{Li}]$ at RPTs for calendar aging cells at multiple temperatures
+### Output of Tuning:
+$$P_\mathrm{cal} = \left[k_{0,\mathrm{SEI}},D_\mathrm{SEI},E_{a,\mathrm{SEI}}\right]^T $$
+
+
+## Step 2
+1. After step 1, we tune the mechanical damage model parameters and lithium plating model parameters together using the cycling aging data at multiple C-rates
+2. The vector of degradation mechanism parameters P to be tuned in step 2 are:
+$$P_\mathrm{cyc} = \left[\beta^-_{\mathrm{LAM},1},\beta^-_{\mathrm{LAM},2},\beta^+_{\mathrm{LAM},1},\beta^+_{\mathrm{LAM},2},m_{\mathrm{LAM}},k_\mathrm{pl}\right] $$
+3. The vector of data to be fitted against model outputs are:
+$$Y_\mathrm{cyc} = \left[C, n_\mathrm{Li},C_n,C_p\right]^T$$
+4. The parameters $P_{cyc}$ were tuned based on the data sets $Y_{cyc}$ extracted from the RPTs of the cycling cells at charge-discharges rates of C/5-C/5%–100%DOD, 1.5C-1.5C-100%DOD, and C/5-1.5C-100%DOD. The reason for tuning the model using the three cells together is to establish the C-rate dependency of the mechanical damage and Li plating models. Using data from a single C-rate would be insufficient to get the C-rate dependency correct in the model. 
+5. Please note that step 2 has to be run after step 1 and we use the values of $P_{cal}$ obtained in step 1 to initialize the model in step 2. 
+
+To perform step 2, please run this [notebook](./step_2_cycling.ipynb)
+
+### Input Data Required:
+eSOH parameters $[x_0,x_{100},y_0,y_{100},C_n,C_p,C,n_{Li}]$ at RPTs for cycling aging cells at multiple C-rates
+### Output of Tuning:
+$$P_\mathrm{cyc} = \left[\beta^-_{\mathrm{LAM},1},\beta^-_{\mathrm{LAM},2},\beta^+_{\mathrm{LAM},1},\beta^+_{\mathrm{LAM},2},m_{\mathrm{LAM}},k_\mathrm{pl}\right] $$
+
+## Step 3
+
+## Step 4

degradation_model/tuning_code/step_1_calendar.ipynb

@@ -0,0 +1,287 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pybamm\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import dfols\n",
+    "from scipy import interpolate\n",
+    "import os, sys\n",
+    "sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(\"__file__\"))))\n",
+    "from batfuns import *\n",
+    "import pickle\n",
+    "plt.rcParams = set_rc_params(plt.rcParams)\n",
+    "\n",
+    "eSOH_DIR = \"../data/esoh/\"\n",
+    "oCV_DIR = \"../data/ocv/\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get parameter values\n",
+    "parameter_values = get_parameter_values()\n",
+    "# define model\n",
+    "spm = pybamm.lithium_ion.SPM(\n",
+    "    {\n",
+    "        \"SEI\": \"ec reaction limited\",\n",
+    "        \"loss of active material\": \"stress-driven\",\n",
+    "        \"lithium plating\": \"irreversible\",\n",
+    "        \"stress-induced diffusion\": \"false\",\n",
+    "    }\n",
+    ")\n",
+    "# spm.print_parameter_info()\n",
+    "param=spm.param\n",
+    "# define length of 1 cycle as 1 day i.e. 24 hours\n",
+    "calendar_time = 24\n",
+    "# initial conditions for ksei, Dsei, Ea,sei\n",
+    "ic= [1.5e-16, 2e-18, 10000]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Tuning"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define the cost function using the vector of data to be fitted\n",
+    "def multi_objective(model, data):\n",
+    "    days = np.floor(data[\"Time\"])\n",
+    "    # Y_cal\n",
+    "    variables = [\"Capacity [A.h]\", \"n_Li\",\"y_0\",\"x_100\"]\n",
+    "    # weights\n",
+    "    weights = [1/5,5,1,1]\n",
+    "    return np.concatenate([\n",
+    "        (np.array(model.loc[days][var]) - np.array(data[var])) * w\n",
+    "        for w,var in zip(weights,variables)\n",
+    "    ]\n",
+    "    )\n",
+    "\n",
+    "# code to simulate lifetime aging\n",
+    "def simulate(x,eps_n_data,eps_p_data,SOC_0,Temp,experiment):\n",
+    "    parameter_values.update(\n",
+    "        {   \n",
+    "            \"Negative electrode LAM min stress [Pa]\": 0,\n",
+    "            \"Negative electrode LAM max stress [Pa]\": 0,\n",
+    "            \"Positive electrode LAM min stress [Pa]\": 0,\n",
+    "            \"Positive electrode LAM max stress [Pa]\": 0,\n",
+    "            \"Positive electrode LAM constant proportional term [s-1]\": 0,\n",
+    "            \"Negative electrode LAM constant proportional term [s-1]\": 0,\n",
+    "            \"Positive electrode LAM constant proportional term 2 [s-1]\": 0,\n",
+    "            \"Negative electrode LAM constant proportional term 2 [s-1]\": 0,\n",
+    "            \"Positive electrode LAM constant exponential term\": 1,\n",
+    "            \"Negative electrode LAM constant exponential term\": 1,\n",
+    "            \"Lithium plating kinetic rate constant [m.s-1]\": 0,\n",
+    "            \"Li plating resistivity [Ohm.m]\": 0,\n",
+    "            \"SEI kinetic rate constant [m.s-1]\": x[0] * ic[0],\n",
+    "            \"EC diffusivity [m2.s-1]\": x[1] * ic[1],\n",
+    "            \"SEI growth activation energy [J.mol-1]\": x[2] * ic[2],\n",
+    "            \"SEI resistivity [Ohm.m]\": 30000.0,\n",
+    "            \"Initial inner SEI thickness [m]\": 0e-09,\n",
+    "            \"Initial outer SEI thickness [m]\": 5e-09,\n",
+    "            \"Negative electrode active material volume fraction\": eps_n_data,\n",
+    "            \"Positive electrode active material volume fraction\": eps_p_data,\n",
+    "            \"Initial temperature [K]\": 273.15+Temp,\n",
+    "            \"Ambient temperature [K]\": 273.15+Temp,\n",
+    "            \"Negative electrode partial molar volume [m3.mol-1]\": 7e-06,\n",
+    "\n",
+    "        },\n",
+    "        check_already_exists=False,\n",
+    "    )\n",
+    "    return cycle_adaptive_simulation_V2(spm, parameter_values, experiment,SOC_0,save_at_cycles=1)\n",
+    "\n",
+    "# function to generate the prediction error\n",
+    "def prediction_error(x):\n",
+    "    try:\n",
+    "        out=[]\n",
+    "        # cell 22: Hot Temp (45C) Calendar Aging @ 100% SOC\n",
+    "        # cell 24: Cold Temp (-5C) Calendar Aging @ 100% SOC\n",
+    "        for cell in [22,24]:\n",
+    "            # load the data for each cell\n",
+    "            cell_no,dfe,N,dfo_0 = load_data_calendar(cell,eSOH_DIR,oCV_DIR)\n",
+    "            # get parameters to intialize the model\n",
+    "            eps_n_data,eps_p_data,SOC_0,Temp = init_exp_calendar(cell_no,dfe,param,parameter_values)\n",
+    "            # define the calendar aging experiment\n",
+    "            experiment = pybamm.Experiment(\n",
+    "                [\n",
+    "                    (\"Rest for \"+f'{calendar_time}'+\" hours\",)\n",
+    "                ]*int(N[-1]+10) ,\n",
+    "                termination=\"50% capacity\",\n",
+    "            #     cccv_handling=\"ode\",\n",
+    "            )\n",
+    "            # simulate the model\n",
+    "            model = simulate(x,eps_n_data,eps_p_data,SOC_0,Temp,experiment)\n",
+    "            # generate the error using cost function\n",
+    "            out_t =   multi_objective(pd.DataFrame(model), dfe)\n",
+    "            # concatenate the error for two cells\n",
+    "            out=np.concatenate([out,out_t])\n",
+    "        print(f\"x={x}, norm={np.linalg.norm(out)}\")\n",
+    "    except Exception as error:\n",
+    "        # error handling\n",
+    "        print(error)\n",
+    "        out=[]\n",
+    "        for cell in [22,24]:\n",
+    "            cell_no,dfe,N,dfo_0 = load_data_calendar(cell,eSOH_DIR,oCV_DIR)\n",
+    "            out_t = np.concatenate([np.array(dfe['Cap'])]*4)\n",
+    "            out=np.concatenate([out,out_t])\n",
+    "        out = 2*np.ones_like(out)\n",
+    "        print(f\"x={x}, norm={np.linalg.norm(out)}\")\n",
+    "    return out"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# initial guess\n",
+    "x0 = np.array([1.0,1.0,1.0])\n",
+    "# lower bound\n",
+    "lower = np.array([1e-3, 1e-3, 1e-3])\n",
+    "# upper bound\n",
+    "upper = np.array([1e+3, 1e+3, 1e+3])\n",
+    "# run tuning\n",
+    "soln_dfols = dfols.solve(prediction_error, x0,bounds=(lower, upper), rhoend=1e-2) #rhoend=1e-4"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# extract solution\n",
+    "x = soln_dfols.x\n",
+    "# multiply back the initial condition to the solution\n",
+    "x1 = np.multiply(ic,x)\n",
+    "# export to dataframe\n",
+    "df_x = pd.DataFrame(columns=['x_0','x_1','x_2','obj'], index=[0])\n",
+    "df_x['x_0'][0]= x1[0]\n",
+    "df_x['x_1'][0]= x1[1]\n",
+    "df_x['x_2'][0]= x1[2]\n",
+    "df_x['obj'][0]= np.linalg.norm(prediction_error(x))\n",
+    "# save results\n",
+    "df_x.to_csv(\"step_1_results.csv\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Simulate Results to Verify"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plotting function\n",
+    "def plotc(all_sumvars_dict1,all_sumvars_dict2,esoh_data,esoh_data_2,leg1=\"sim1\",leg2=\"sim2\"):\n",
+    "    esoh_vars = [\"x_100\", \"y_0\", \"C_n\", \"C_p\", \"Capacity [A.h]\", \"Loss of lithium inventory [%]\"]\n",
+    "    fig, axes = plt.subplots(3,2,figsize=(7,7))\n",
+    "    for k, name in enumerate(esoh_vars):\n",
+    "        ax = axes.flat[k]\n",
+    "        ax.plot(all_sumvars_dict1[\"Cycle number\"],all_sumvars_dict1[name],\"r^\")\n",
+    "        ax.plot(esoh_data[\"N\"],esoh_data[name],\"kx\")\n",
+    "        ax.plot(all_sumvars_dict2[\"Cycle number\"],all_sumvars_dict2[name],\"bv\")\n",
+    "        ax.plot(esoh_data_2[\"N\"],esoh_data_2[name],\"kx\")\n",
+    "        ax.set_title(split_long_string(name))\n",
+    "        if k ==2 or k==3:\n",
+    "            ax.set_ylim([5,6.2])\n",
+    "        if k>3:\n",
+    "            ax.set_xlabel(\"Days\")\n",
+    "    fig.legend([leg1, \"Data\", leg2 ], \n",
+    "           loc=\"lower center\",bbox_to_anchor=[0.5,-0.05], ncol=1, fontsize=11)\n",
+    "    fig.tight_layout()\n",
+    "    return fig"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol = {}\n",
+    "esoh_data = {}\n",
+    "i = 0\n",
+    "# cell 22: Hot Temp (45C) Calendar Aging @ 100% SOC\n",
+    "# cell 24: Cold Temp (-5C) Calendar Aging @ 100% SOC\n",
+    "for cell in [22,24]:\n",
+    "    # load the data for each cell\n",
+    "    cell_no,dfe,N,dfo_0 = load_data_calendar(cell,eSOH_DIR,oCV_DIR)\n",
+    "    esoh_data[i]=dfe\n",
+    "    # get parameters to intialize the model\n",
+    "    eps_n_data,eps_p_data,SOC_0,Temp = init_exp_calendar(cell_no,dfe,param,parameter_values)\n",
+    "    # define the calendar aging experiment\n",
+    "    experiment = pybamm.Experiment(\n",
+    "        [\n",
+    "            (\"Rest for \"+f'{calendar_time}'+\" hours\",)\n",
+    "        ]*int(N[-1]+10) ,\n",
+    "        termination=\"50% capacity\",\n",
+    "    #     cccv_handling=\"ode\",\n",
+    "    )\n",
+    "    # simulate the model\n",
+    "    sol[i] = simulate(x,eps_n_data,eps_p_data,SOC_0,Temp,experiment)\n",
+    "    i+=1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot tuning results\n",
+    "# plot hot and cold data and simulations together\n",
+    "plotc(sol[0],sol[1],esoh_data[0],esoh_data[1],leg1=\"hot\",leg2=\"cold\");"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.9.7 ('windows-dev': venv)",
+   "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.9.13"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "fe6e1a396f7757d05554318f98bb2d7e7d3785df9fdbcce0707f057fad4349a9"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}