{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Full-cell OCV\n",
    "\n",
    "In this example, we use synthetic full-cell GITT data and known half-cell OCP parameters to determine the stoichiometry windows that give the correct full-cell OCV, a process often referred to as \"electrode balancing\".\n",
    "\n",
    "Before running this example, make sure to run the script `true_parameters/generate_data.py` to generate the synthetic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import ionworkspipeline as iwp\n",
    "import iwutil\n",
    "\n",
    "from true_parameters.parameters import get_msmr_parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cell balancing theory \n",
    "\n",
    "Given the open-circuit potentials for each electrode, we combine them to get the full-cell open-circuit voltage\n",
    "$$U_{cell}(z) = U_{p}(\\theta_{p}) - U_{n}(\\theta_{n}),$$\n",
    "where $z$ is the full-cell state of charge, and $\\theta_{n}$ and $\\theta_{p}$ are the electrode-level state of charge of the negative and positive electrodes respectively (usually called \"stoichiometry\" or \"lithiation\" to avoid confusion with cell-level state of charge). To evaluate the OCV, we need to convert the stoichiometry of each electrode to the full-cell state of charge. This is done using the following formula:\n",
    "$$\n",
    "    z = \\frac{\\theta_{n} - \\theta_{n}^\\mathrm{0}}{\\theta_{n}^\\mathrm{100} - \\theta_{n}^\\mathrm{0}} = \\frac{\\theta_{p} - \\theta_{p}^\\mathrm{0}}{\\theta_{p}^\\mathrm{100} - \\theta_{p}^\\mathrm{0}},\n",
    "$$\n",
    "where $\\theta_{n}^\\mathrm{0}$, $\\theta_{n}^\\mathrm{100}$, $\\theta_{p}^\\mathrm{0}$, $\\theta_{p}^\\mathrm{100}$ are the stoichiometries of the negative and positive electrodes at 0\\% and 100\\% state of charge respectively. Hence the challenge is to find the four values $\\theta_{n}^\\mathrm{0}$, $\\theta_{n}^\\mathrm{100}$, $\\theta_{p}^\\mathrm{0}$, $\\theta_{p}^\\mathrm{100}$, that best fit the full-cell open-circuit voltage data.\n",
    "In practice, full-cell OCV data is is given in terms of capacity instead of stoichiometry. Therefore, we reformulate the full-cell open-circuit voltage as a function of capacity:\n",
    "$$U_{cell}(q) = U_{p}(q_{p}) - U_{n}(q_{n}),$$\n",
    "where $q = Qz$, $q_{n} = \\theta_{n}Q_{n}$, $q_{p} = \\theta_{p}Q_{p}$, and $Q$ is the total capacity of the cell and $Q_{n}$ and $Q_{p}$ are the total capacities of the individual electrodes. We can then convert between the full-cell capacity and the electrode capacities using the following formula:\n",
    "$$\n",
    "    q/Q = \\frac{q_{n} - q_{n}^\\mathrm{0}}{q_{n}^\\mathrm{100} - q_{n}^\\mathrm{0}} = \\frac{q_{p} - q_{p}^\\mathrm{0}}{q_{p}^\\mathrm{100} - q_{p}^\\mathrm{0}}.\n",
    "$$\n",
    "We find the values of $q_{n}^\\mathrm{0}$, $q_{n}^\\mathrm{100}$, $q_{p}^\\mathrm{0}$, $q_{p}^\\mathrm{100}$ that best fit the data.\n",
    "\n",
    "In this example, we fit the \"lower excess capacity\" and \"upper excess capacity\" instead of the stoichiometries at 0\\% and 100\\% state of charge. We use these names because the data typically does not reach the true 0\\% and 100\\% state of charge, due to kinetic limitations and/or electrolyte stability. The lower excess capacity is equal to $q_{n}^\\mathrm{0}$ in the negative electrode, and $q_{p}^\\mathrm{100}$ in the positive electrode. The upper excess capacity is equal to $Q_{n} - q_{n}^\\mathrm{100}$ in the negative electrode, and $Q_{p} - q_{p}^\\mathrm{0}$ in the positive electrode. $Q$ is given by the actual total capacity observed in the experimental data (called \"useable capacity\" in the example)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data\n",
    "\n",
    "We load the synthetic GITT data and the extract the OCP from the relaxed voltages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gitt_data = iwutil.read_df(Path(\"synthetic_data\") / \"full_cell\" / \"gitt.csv\")\n",
    "ocp_data = iwp.data_fits.preprocess.pulse_data_to_ocp(gitt_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the known parameters\n",
    "\n",
    "In this example, we assume we already know the MSMR OCP parameters for the negative and positive electrodes. In practice, these parameters can be fitted using the half-cell OCV workflow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msmr_params_n = get_msmr_parameters(\"negative\")\n",
    "msmr_params_p = get_msmr_parameters(\"positive\")\n",
    "known_params = {\n",
    "    **msmr_params_n,\n",
    "    **msmr_params_p,\n",
    "    \"Ambient temperature [K]\": 298.15,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the parameters to fit\n",
    "\n",
    "Now we set up our initial guesses for the parameters. As in the previous example, we create a dictionary of `Parameter` objects.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_use = ocp_data[\"Capacity [A.h]\"].max()\n",
    "Q_n_lowex = iwp.Parameter(\n",
    "    \"Q_n_lowex\", initial_value=0.01 * Q_use, bounds=(0, 0.2 * Q_use)\n",
    ")\n",
    "Q_p_lowex = iwp.Parameter(\"Q_p_lowex\", initial_value=0.1 * Q_use, bounds=(0, Q_use))\n",
    "Q_n_uppex = iwp.Parameter(\n",
    "    \"Q_n_uppex\", initial_value=0.1 * Q_use, bounds=(0, 0.2 * Q_use)\n",
    ")\n",
    "Q_p_uppex = iwp.Parameter(\"Q_p_uppex\", initial_value=0.5 * Q_use, bounds=(0, Q_use))\n",
    "parameters = {\n",
    "    \"Negative electrode lower excess capacity [A.h]\": Q_n_lowex,\n",
    "    \"Positive electrode lower excess capacity [A.h]\": Q_p_lowex,\n",
    "    \"Negative electrode upper excess capacity [A.h]\": Q_n_uppex,\n",
    "    \"Positive electrode upper excess capacity [A.h]\": Q_p_uppex,\n",
    "    \"Negative electrode capacity [A.h]\": Q_n_lowex + Q_use + Q_n_uppex,\n",
    "    \"Positive electrode capacity [A.h]\": Q_p_lowex + Q_use + Q_p_uppex,\n",
    "    \"Usable capacity [A.h]\": Q_use,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the workflow\n",
    "\n",
    "Now we are ready to run the workflow. To run the workflow with the default settings, we can use the `fit_plot_save` function. This function runs the fit, saves the results (if a `save_dir` keyword argument is provided), and plots the results. All of the workflows have a common interface, where you typically need to provide the data, the parameters to be fitted, and a dictionary of known parameters.\n",
    "\n",
    "For the full-cell OCV workflow, we need to specify the voltage range of each electrode to use in the fit. We can do this by passing a dictionary of keyword arguments, `objective_kwargs`, to be used in the objective. The limits are used to evaluate the half-cell OCPs in the fit, and should therefore span a wider half-cell voltage range than we expect to see in the data.\n",
    "\n",
    "The workflow returns a `Results` object, which contains the fitted parameters, and a dictionary of figures with the keys \"trace\" and the name of the objective function (in this case \"MSMRFullCell\"). The trace figure shows the evolution of the cost function and parameter values during the fit. The MSMRFullCell figure shows the fitted OCV and dU/dQ curves. In the plot, we can see how the half-cell OCPs combine to give the full-cell OCV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options = {\n",
    "    # these are the voltage ranges over which to evaluate the half-cell OCP curves\n",
    "    \"negative voltage limits\": (0, 2),\n",
    "    \"positive voltage limits\": (2.5, 5),\n",
    "}\n",
    "objective_kwargs = {\n",
    "    \"options\": options,\n",
    "}\n",
    "res, figs = iwp.workflows.full_cell_ocv.fit_plot_save(\n",
    "    ocp_data,\n",
    "    parameters,\n",
    "    known_params,\n",
    "    objective_kwargs=objective_kwargs,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take a look at the fitted parameters by accessing the `parameter_values` attribute of the `Results` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.parameter_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's going on? A lower-level interface\n",
    "\n",
    "The workflows are designed to be simple to use and cover most use cases. They can be customized by passing additional keyword arguments. For example, we can pass a dictionary `datafit_kwargs` to specify the cost function or optimizer to be used in the fit. Much more customization is possible by accessing the underlying `DataFit` and `Optimizer` classes directly. Here we will go through the steps of running the fit manually.\n",
    "\n",
    "We can use the same `parameters` and `known_params` as dicitonaries as before, but we need to create the `Objective` and `DataFit` objects manually.\n",
    "\n",
    "Let's start by creating the objective function. We use the `MSMRFullCell` objective function and pass in the data and voltage limits.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = iwp.objectives.MSMRFullCell(\n",
    "    ocp_data,\n",
    "    options={\n",
    "        \"negative voltage limits\": (0, 2),\n",
    "        \"positive voltage limits\": (2.5, 5),\n",
    "    },\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we've created the objective, we can create the `DataFit` object. The `DataFit` object takes in the objective and the parameters to fit. We can also customize the fitting processing by passing in a cost function and selecting an optimizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cost = iwp.costs.Difference()\n",
    "optimizer = iwp.optimizers.ScipyLeastSquares(verbose=1)\n",
    "datafit = iwp.DataFit(objective, parameters=parameters, cost=cost, optimizer=optimizer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we run the fit and plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = datafit.run(known_params)\n",
    "_ = datafit.plot_trace()\n",
    "_ = datafit.plot_fit_results()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}