proof of an implicit solver

diode-model-schematic.png

explicit-pv-model.ipynb

PV-system-implicit-solver.ipynb

PV-system-MPP-implicit-solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example of implicit approach to solving for max power of a PV system\n",
    "\n",
    "This example will demonstrate the following:\n",
    "\n",
    "1. the design pattern for setting up a numerical solver that takes a callback function\n",
    "2. specifically how to write a callback for the unknowns and constraints to find the max power of a PV system\n",
    "3. using a solver with the callback to find a solution\n",
    "\n",
    "This next section just sets the notebook up to show plots inline and imports NumPy to do vector math and the numerical solvers we use from the SciPy `optimize` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "# NOTE: fsolve is the same as root() set to 'hybr'\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Design patterns for numerical solvers of implicit systems of equations\n",
    "\n",
    "Setting up a numerical solver is a little counter intuitive, because the solver takes a callback as an argument and an initial guess of the optimal solution. In this example I am going to use a root-finding solver, so the optimum solution is the set of callback arguments that return zero.\n",
    "\n",
    "### Implicit systems of equations\n",
    "\n",
    "I think of an implicit system of equations as non-linear. I think of a linear system of equations as something like ...\n",
    "\n",
    "$$Y_i = A_{ij} X_j$$\n",
    "\n",
    "... where you can invert $A_{ij}$ to get $X_j$. In an implicit system, you can't write an _explicit_ expression for any individual $y$ as a function of $x$ because $x$ is also a function of $y$. An example of an implicit equation is the single diode model of a solar cell in which the cell current, $I$, is on both sides:\n",
    "\n",
    "$$I = I_L - I_0 \\left( \\exp \\left( \\frac{V + I R_s}{\\gamma V_t} \\right) - 1 \\right) - \\frac{V + I R_s}{R_{sh}}$$\n",
    "\n",
    "This equation can't be easily written explicitly for $I$. You could try to use a Taylor series approximation to solve for $I$ but it depends on where you expand the approximation, and how many terms from the series you select. Another approach would be to pick a voltage to solve for and then iteratively guess a value for $I$ until the both sides of the equation were equal. This is what a numerical solver does, and it uses a couple of tricks to choose iteratively better guess to find the solution in as few iterations as possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CONSTANTS\n",
    "IL = 7.0  # [A] photogenerated current\n",
    "I0 = 1.0e-6  # [A] dark current\n",
    "RSH = 20.0  # [ohms] shunt resistance\n",
    "RS = 0.001  # [ohms] series resistance\n",
    "GAMMA = 1.23  # diode ideality\n",
    "VT = 0.026  # [V] thermal voltage at 300[K]\n",
    "VBYPASS = -0.5  # bypass diode trigger voltage [V]\n",
    "NSERIES_CELLS = 72  # number of series cells per module\n",
    "NMODS = 8  # number of modules per string\n",
    "NSTRINGS = 3  # number of strings in system\n",
    "NBYPASS = 3  # number of bypass diodes per module, assumes cells divided evenly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Single Diode Model (SDM)\n",
    "\n",
    "The single diode model is an electrical analog of a solar cell.\n",
    "\n",
    "![single diode model schematic](./diode-model-schematic.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sdm_residual_and_der(il, i0, rsh, rs, gamma, vt, vcells, istrings):\n",
    "    \"\"\"\n",
    "    Residual for single diode model (SDM) and derivative.\n",
    "\n",
    "    Args:\n",
    "        il (numeric): photogenerated current [A]\n",
    "        i0 (numeric): dark current [A]\n",
    "        rsh (numeric): shunt resistance [Ohms]\n",
    "        rs (numeric): series resistance [Ohms]\n",
    "        gamma (numeric): diode ideality factor\n",
    "        vt (numeric): thermal voltage [V]\n",
    "        vcells (numeric): cell voltages [V]\n",
    "        istrings (numeric): string currents [A]\n",
    "\n",
    "    Returns:\n",
    "        residuals (numeric): difference between given and calculated ``istring`` [A]\n",
    "        di_dv (numeric): derivative of current versus voltage [A/V]\n",
    "    \"\"\"\n",
    "    vdiode = vcells + istrings*rs  # voltage drop across diode in SDM [V]\n",
    "    nvt = gamma * vt  # product of the ideality factor and thermal voltage [V]\n",
    "    # temporary calculations to simplify the expressions below\n",
    "    a = np.exp(vdiode / nvt)\n",
    "    b = i0*a/nvt + 1/rsh\n",
    "    di_dv = -b / (1.0 + rs*b)  # the derivative of current versus voltage\n",
    "    # calculate the difference between the given and calculated current \n",
    "    residuals = il - i0*(a - 1.0) - vdiode/rsh - istrings\n",
    "    return residuals, di_dv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### unknowns and constraints\n",
    "\n",
    "This trivial example demonstrates finding the current of a single solar cell for a given voltage. If we want to know how an entire PV system performs then there are many unknowns and therefore there will also be many equations. I call these equations **constraints**, and there must be as many constraints as there are unknowns for the system to have a unique solution. Although there's no guaranty that there's a solution at all. Assuming we have the correct set of constraints to find a solution, I call this a square or closed system. If the system is not closed or square then you could have either an over or under constrained problem which is also very interesting, but isn't the topic of this example.\n",
    "\n",
    "## Solving for max power of a PV system\n",
    "\n",
    "A PV system is an array of solar cells that are organized serially into panels and strings and then the strings of panels are arranged in parallel. This arrangement poses some additional constraints, due to Ohm's and Kirchhoff's laws, in addition to the single diode model introduced above.\n",
    "\n",
    "1. The cells in series in a panel must all have the same current.\n",
    "2. The panels in series in a string must also have the same current.\n",
    "3. The strings in parallel in the array must have the same voltage.\n",
    "\n",
    "Before we get to those constraints, let's focus on the desired solution. We want to know what voltage to operate the PV system at in order to get the maximum power. We can write this constraint by realizing that the maximum power occurs where the slope of power as a function of voltage is zero, as indicated by the flat green line in the figure below.\n",
    "\n",
    "![max power point](./solar-convex-optimum.png)\n",
    "\n",
    "Since we are using a root-finding solver, this constraint is perfect because roots are where the function returns zero. So our callback should return the derivative $\\frac{dP_{sys}}{dV_{sys}}$ which will be zero at the max power. The system power, $P_{sys}$ and voltage, $V_{sys}$, are both unknowns, so our system is not closed and we need more constraints. So we can rewrite the derivate in terms of the system current:\n",
    "\n",
    "$$\\frac{dP_{sys}}{dV_{sys}} = \\frac{d\\left(I_{sys}V_{sys}\\right)}{dV_{sys}} = I_{sys} + V_{sys} \\frac{dI_{sys}}{dV_{sys}}$$\n",
    "\n",
    "Now we have three unknowns, and still only one constraint. So let's use what we know about the single diode model and how a PV system is arranged to generate the rest of the constraints. First let's use Kirchhoff's law to expand $I_{sys}$ which is the sum of the string currents, since they're in parallel.\n",
    "\n",
    "$$\\frac{dP_{sys}}{dV_{sys}} = \\sum{I_{str}} + V_{sys} \\frac{d \\sum{I_{str}}}{dV_{sys}}$$\n",
    "\n",
    "Then let's take the trivial case of two strings each with two 60-cell panels (or _modules_ as they are usually called), and let's just look at the derivative of $\\frac{d \\sum{I_{str}}}{dV_{sys}}$.\n",
    "\n",
    "$$\\frac{d \\sum{I_{str}}}{dV_{sys}} = \\frac{d \\left(I_{str1} + I_{str2}\\right)}{dV_{sys}}\n",
    "  = \\frac{d I_{str1}}{dV_{sys}} + \\frac{d I_{str2}}{dV_{sys}}$$\n",
    "\n",
    "Since these are definite derivates not partials we can write their reciprocal as $\\frac{dI}{dV}=\\frac{1}{\\frac{dV}{dI}}$.\n",
    "\n",
    "$$\\frac{d \\sum{I_{str}}}{dV_{sys}} = \\frac{1}{\\frac{dV_{sys}}{d I_{str1}}} + \\frac{1}{\\frac{dV_{sys}}{d I_{str2}}}$$\n",
    "\n",
    "Now let's use Ohm's law to expand $V_{sys}$ which is the sum of the cell voltages for each string since they're in series.\n",
    "\n",
    "$$\\frac{d \\sum{I_{str}}}{dV_{sys}} = \\frac{1}{\\frac{d \\sum{V_{1j}}}{d I_{str1}}} + \\frac{1}{\\frac{d \\sum{V_{2j}}}{d I_{str2}}}\n",
    "  = \\frac{1}{\\sum{\\frac{d V_{1j}}{d I_{str1}}}} + \\frac{1}{\\sum{\\frac{d V_{2j}}{d I_{str2}}}}$$\n",
    "\n",
    "Recall that the cell current $I$ will be the same for all cells in series in the string, so all the cells in the 1<sup>st</sup> string have $I_{str1}$, $I_{str2}$ in the 2<sup>nd</sup> string, and so on for each string $i$. Therefore we can rewrite the constraint on max power as follows:\n",
    "\n",
    "$$\\frac{dP_{sys}}{dV_{sys}} = \\sum{I_{str}} + V_{sys} \\frac{1}{\\sum{\\frac{d V_{ij}}{d I_{str,i}}}}\n",
    " = \\sum{I_{str}} + V_{sys} \\sum{ \\frac{1}{\\sum{ \\frac{1}{ \\frac{d I_{str,i}}{d V_{ij}} } }} }$$\n",
    "\n",
    "It's easier to calculate $\\frac{dI}{dV}$ so $\\frac{dV}{dI}$ was rewritten as the reciprocal in the expression above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dp_dv(istrings, nstrings, di_dv, vsys):\n",
    "    \"\"\"\n",
    "    Slope of PV system power-voltage curve\n",
    "\n",
    "    Args:\n",
    "        istrings (numeric): string currents [A]\n",
    "        nstrings (int): number of strings must be ``len(istrings)``\n",
    "        di_dv (numeric): derivatives of current versus cell voltages,\n",
    "            must have shape ``[nseries_cells * nmods, nstrings]``\n",
    "        vsys (float): system voltage\n",
    "\n",
    "    Returns:\n",
    "        dp_dv (float): the slope of the power-voltage curve\n",
    "\n",
    "\n",
    "    Use this as the residual for max power point (MPP) since dP/dV = 0 at MPP.\n",
    "    \"\"\"\n",
    "    isys = np.sum(istrings)\n",
    "    dv_di = [np.sum(1/di_dv[:, string]) for string in range(nstrings)]\n",
    "    return isys + vsys * np.sum(1 / np.array(dv_di))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have all of the cell voltages, $V_{ij}$, the string currents $I_{str,i}$, and the system voltage $V_{sys}$ as unknowns. We can use the single diode model equation as the constraint for the cell voltages\n",
    "\n",
    "$$F_{ij} = I_L - I_0 \\left( \\exp \\left( \\frac{V_{ij} + I_j R_s}{\\gamma V_t} \\right) - 1 \\right) - \\frac{V_{ij} + I_j R_s}{R_{sh}} - I_j = 0$$\n",
    "\n",
    "... and we can use the sum of the cell voltages as the constraint for he string currents\n",
    "\n",
    "$$F_i = V_{sys} - \\sum{V_{ij}}$$\n",
    "\n",
    "... and our system of equations is closed!\n",
    "\n",
    "#### writing the callback\n",
    "Now we can write our callback. This is a function that gets called with each updated guess everytime the solver iterates. The form of the callback function is specified by the solver, but it almost always has the form `callback(x, *args)` where `x` is a vector of the guesses and `args` is a tuple of a extra arguments that are constant as far as the solver is concerned. The solver we are using in this example is called [`fsolve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html), it's part of the SciPy optimize package, and it's a wrapper around the MINPACK solver called HYBRJ which is a modified version of the Powell method. Writing the callback function is basically book-keeping, since we just need to account for the order of the inputs `args`, outputs `x`, and the constraints which are evaluated and returned as residuals.\n",
    "\n",
    "##### extra arguments\n",
    "The first part of the callback expands the extra arguments and gives them names that are more convenient. Since `args` is a tuple their order is important. These arguments are ones that we want to change between different scenarios, but are considered constant to to the solver. For example we are going to solve a PV system with 3 strings, 8 modules per string, 72 cells per module, and 3 bypass diodes per module, which correspond to the 3<sup>rd</sup>, 2<sup>nd</sup>, 1<sup>st</sup>, and 4<sup>th</sup> values in `args`. To solve a different PV array configuration we would just change values in `args`.\n",
    "\n",
    "Another way to think of `args` is as the independent variables, and the guesses, `x`, are the dependent variables. This may seem counterintuitive, since the solver passes both the guesses and extra arguments to the callback. Normally function returns the dependent variables given the independent variabls. This is why an implicit solution is sometimes called a reverse problem. The other independent variables that we're going to put into `args` are the photogenerated current, dark current, shunt resistance, series resistance, diode ideality factor, thermal voltage, and the bypass trigger voltage.\n",
    "\n",
    "##### dependent variables\n",
    "The next part of the callback expands the first argument, `x`, which are the unknown dependent variables that we are solving for. Each iteration the solver makes a guess of what it thinks the solution should be, passes `x` to the callback, and checks if the residuals are close to zero. The order of the dependent variables is important. We unpack the cell voltages first and then reshape them so that each column is a string, this makes it easier to iterate over the strings and aggregate the voltages. Then we unpack the submodule voltages, string currents, and finally the system voltage, which is the value we're really interested in remember.\n",
    "\n",
    "##### constraints\n",
    "After all of the inputs `args` and outputs `x` are unpacked and given convenient names and reshaped, we evaluate the constriants and assemble the residuals to return to the solver. The order of the residuals may influence the solver convergence rate, but this is a topic for a different example. Our problem is actually very sparse, meaning most of the constriants only depend on a few dependent variables, and so the best way to organize the residuals so that the Jacobian (derivative of constraints with respect to `x`) is diagonal, which we can do by evaluating the constraints in the same order as `x`.\n",
    "\n",
    "The first constraint we evaluate are the single diode model residuals given by `sdm_residual_and_der` which we wrote above. This is a vector function which we pass in the single diode model parameters, the cell voltages, and the string currents, and returns an array of residuals with the same shape as the cell voltages. For convenience, we've also written this constraint to return the derivative $\\frac{dI}{dV}$ so we can use again later.\n",
    "\n",
    "###### bypass diodes\n",
    "The next constraint is actually one we haven't talked about yet. It checks the submodule voltages to see if they are low enough to trigger the bypass diode that is in parallel to the submodule and limits the amount of power lossed when a cell is in reverse bias.\n",
    "\n",
    "The last two residuals are given by a constriant on the system voltage, which should be the sum of the cell voltages for each string, and the slope of the power verses voltage curve, which will be zero at the max power point. We use the function written above called `dp_dv` to evaluate that last constraint, append it to the end of the residuals, and then finally return all of the residuals to the solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mpp(x, *args):\n",
    "    \"\"\"\n",
    "    Callback function for `fsolve` to find max power point of a PV system.\n",
    "    \n",
    "    Args:\n",
    "        x (numeric): dependent variables\n",
    "        args (tuple): independent variables and any extra arguments required\n",
    "            by callback, considered constant by solver\n",
    "\n",
    "    Returns:\n",
    "        residuals (numeric): the constraints evaluated given ``x`` and ``args``\n",
    "\n",
    "    The dependent variables, ``x``,  must be in the following order:\n",
    "    cell voltages, submodule voltages, string currents, and system voltage.\n",
    "\n",
    "    The independent variables, ``args``, must be in the following order:\n",
    "    number of series cells per module, number of modules per string, number of\n",
    "    strings, number of bypass diodes per module, and single diode parameters.\n",
    "\n",
    "    The single diode parameters must be in the following order:\n",
    "    photogenerated current, dark current, shunt resistance, series resistance,\n",
    "    diode ideality factor, thermal voltage, and bypass diode trigger voltage.\n",
    "    \"\"\"\n",
    "    # expand extra arguments\n",
    "    nseries_cells = args[0]  # number of series cells per module\n",
    "    nmods = args[1]  # number of modules per string\n",
    "    nstrings = args[2]  # number of strings\n",
    "    # number of bypyass diodes per module, assumes cells divided evenly\n",
    "    nbypass = args[3]\n",
    "\n",
    "    # single diode parameters and module bypass diode trigger voltage\n",
    "    il, i0, rsh, rs, gamma, vt, vbypass = args[4:]\n",
    "\n",
    "    # expand the current guesses for cell voltages, submodule voltages,\n",
    "    # string currents, and system voltage\n",
    "\n",
    "    # cell voltages\n",
    "    idx0 = nseries_cells * nmods * nstrings\n",
    "    # reshape vcells as 2-D so 2nd dimension are strings\n",
    "    vcells = x[:idx0].reshape(-1, nstrings)\n",
    "\n",
    "    # submodule voltages\n",
    "    idx1 = idx0 + nbypass * nmods * nstrings\n",
    "    # reshape vsubs as 2-D so 2nd dimension are strings\n",
    "    vsubs = x[idx0:idx1].reshape(-1, nstrings)\n",
    "\n",
    "    # strings currents\n",
    "    idx2 = idx1 + nstrings\n",
    "    istrings = x[idx1:idx2]\n",
    "\n",
    "    # system voltage\n",
    "    vsys = x[-1]\n",
    "\n",
    "    # residuals\n",
    "    residuals = np.zeros(x.shape)\n",
    "\n",
    "    # evalulate the constraints using the current guesses, and organize the\n",
    "    # residuals to return to the solver\n",
    "\n",
    "    # cells constraints\n",
    "    sdm_residual, di_dv = sdm_residual_and_der(il, i0, rsh, rs, gamma, vt, vcells, istrings)\n",
    "    residuals[:idx0] = np.ravel(sdm_residual)\n",
    "\n",
    "    # bypass diode constraint\n",
    "    ncells_bypass = nseries_cells / nbypass  # assume cells are divided evenly\n",
    "    vsubmodule = np.zeros(vsubs.shape)  # allocate space for submodule voltages\n",
    "    # loop over strings and submodules to add up their the cell voltages\n",
    "    for string in range(nstrings):\n",
    "        for submodule in range(nbypass * nmods):\n",
    "            jdx = int(submodule*ncells_bypass)\n",
    "            vsubmodule[submodule, string] = np.sum(vcells[jdx:int(jdx+ncells_bypass), string])\n",
    "    # activate the bypass diode if the submodule voltage exceeds the trigger voltage\n",
    "    vsubmodule[vsubmodule < vbypass] = vbypass\n",
    "    residuals[idx0:idx1] = np.ravel(vsubmodule - vsubs)\n",
    "\n",
    "    # string constraints\n",
    "    vstrings = [np.sum(vsubs[:, string]) for string in range(nstrings)]\n",
    "    residuals[idx1:idx2] = vsys - np.array(vstrings)\n",
    "\n",
    "    # system constraint\n",
    "    residuals[-1] = dp_dv(istrings, nstrings, di_dv, vsys)\n",
    "\n",
    "    return residuals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve it\n",
    "\n",
    "Finally we call the solver with the callback and extra arguments, `args`, and an initial guess to get it started. The initial guess can influence the solver convergence rate, so try to pick something reasonable. We pick some abitrary voltages and currents that should be in the the range of the cell IV curves. I also add some _noise_ to the photogenereated current so that every cell is different to make the solution interesting. Otherwise the problem is trivial. Then we look at the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'nfev': 1826,\n",
       "  'fjac': array([[-9.99626699e-01,  8.39381518e-04,  1.27058288e-04, ...,\n",
       "          -2.87685630e-09, -2.87262033e-09, -4.46879831e-04],\n",
       "         [-7.76128467e-05, -9.99332760e-01, -4.11900480e-04, ...,\n",
       "           9.15969279e-10,  9.25541462e-10, -4.43889477e-04],\n",
       "         [ 6.57934909e-05,  1.06368279e-04, -9.99360978e-01, ...,\n",
       "          -6.78712063e-10, -6.89165936e-10, -4.53175677e-04],\n",
       "         ...,\n",
       "         [ 1.16852415e-04,  1.25416214e-03,  7.99807466e-04, ...,\n",
       "          -1.51937321e-01, -1.16291962e-02, -6.90958484e-02],\n",
       "         [-8.32766976e-06,  8.55342780e-05,  9.31119629e-04, ...,\n",
       "           8.18983682e-02, -1.62966369e-01,  2.05052489e-03],\n",
       "         [ 1.05967505e-03,  1.09058168e-03,  1.04506858e-03, ...,\n",
       "           6.40370881e-02,  6.40365172e-02, -8.15358915e-01]]),\n",
       "  'r': array([ 2.32609894e+02, -1.61211926e-01, -7.99751864e-02, ...,\n",
       "         -7.47145942e-01,  2.68675465e-03,  8.64189840e-01]),\n",
       "  'qtf': array([ 3.19480703e-10,  1.01757301e-10, -1.26874569e-10, ...,\n",
       "          3.15419133e-11, -4.16989199e-11, -2.99086638e-11]),\n",
       "  'fvec': array([ 1.64044778e-10,  1.32882683e-09, -1.35626887e-09, ...,\n",
       "          2.84217094e-14,  2.84217094e-14, -7.67386155e-13])},\n",
       " 1,\n",
       " 'The solution converged.')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# make an initial guess\n",
    "varb = 0.5  # [V] some voltage - use voc_est with Rsh=np.Inf & Rs=0\n",
    "iarb = 6.0  # [A] some current\n",
    "v0 = [varb] * (NSTRINGS*NMODS*NSERIES_CELLS)  # [V] cell voltages guesses\n",
    "vsub0 = [varb * NSERIES_CELLS/NBYPASS] * (NSTRINGS*NMODS*NBYPASS)  # [V] submodule voltages guesses\n",
    "istr0 = [iarb] * NSTRINGS  # [A] string currents guesses\n",
    "vsys0 = varb*NSERIES_CELLS*NMODS\n",
    "X0 = np.array(v0 + vsub0 + istr0 + [vsys0])  # inital guess\n",
    "\n",
    "# add some noise to the light current\n",
    "noise = np.random.randn(NMODS*NSERIES_CELLS, NSTRINGS) / 1000\n",
    "\n",
    "# solve it\n",
    "ARGS = (NSERIES_CELLS, NMODS, NSTRINGS, NBYPASS, IL+noise, I0, RSH, RS, GAMMA, VT, VBYPASS)\n",
    "full_output = fsolve(mpp, x0=X0, args=ARGS, full_output=True)  # if full output then results is a tuple\n",
    "results = full_output[0]\n",
    "full_output[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.41360451, 0.41356066, 0.41368987],\n",
       "       [0.41364739, 0.41358087, 0.41345164],\n",
       "       [0.4135863 , 0.413635  , 0.4135964 ],\n",
       "       ...,\n",
       "       [0.41363315, 0.4136311 , 0.41374838],\n",
       "       [0.41352228, 0.41371391, 0.41358415],\n",
       "       [0.41353419, 0.41364702, 0.41353581]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx0 = NSTRINGS*NMODS*NSERIES_CELLS\n",
    "voltages = results[:idx0].reshape(-1, NSTRINGS)  # cell voltages\n",
    "voltages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9.9265523 , 9.92619247, 9.92649578],\n",
       "       [9.92723563, 9.92716908, 9.92633166],\n",
       "       [9.92611235, 9.92631447, 9.92651405],\n",
       "       [9.9261438 , 9.92611915, 9.92670042],\n",
       "       [9.92687155, 9.92674808, 9.92725245],\n",
       "       [9.92659276, 9.92707929, 9.92682993],\n",
       "       [9.92681792, 9.92707927, 9.92662459],\n",
       "       [9.92680918, 9.92625955, 9.9264328 ],\n",
       "       [9.92636344, 9.92631101, 9.9266266 ],\n",
       "       [9.92683002, 9.92619413, 9.92703828],\n",
       "       [9.92623212, 9.92761628, 9.92655679],\n",
       "       [9.92647733, 9.92712028, 9.92674466],\n",
       "       [9.92730756, 9.92646369, 9.92658204],\n",
       "       [9.92674556, 9.92679508, 9.92659195],\n",
       "       [9.92643695, 9.92599012, 9.92691555],\n",
       "       [9.92676782, 9.92695999, 9.92675666],\n",
       "       [9.927116  , 9.92659797, 9.92636442],\n",
       "       [9.92637441, 9.92723361, 9.92673005],\n",
       "       [9.92701461, 9.92676741, 9.92657251],\n",
       "       [9.92705106, 9.92667462, 9.92674367],\n",
       "       [9.9267415 , 9.92690679, 9.92661515],\n",
       "       [9.92659617, 9.92658591, 9.92696929],\n",
       "       [9.92696229, 9.92646683, 9.92632834],\n",
       "       [9.92631406, 9.92682132, 9.92714877]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idx1 = idx0 + NSTRINGS*NMODS*NBYPASS\n",
    "submodules = results[idx0:idx1].reshape(-1, NSTRINGS)\n",
    "submodules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6.47225027, 6.47225118, 6.47224119])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "strings = results[idx1:(idx1+NSTRINGS)]\n",
    "strings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "238.2404664032513"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = results[-1]\n",
    "system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "19.416742651704602"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isys = sum(strings)\n",
    "isys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4625.853825374007"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pmp = system * isys\n",
    "pmp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time how long it takes to run this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.28 s ± 117 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit fsolve(mpp, x0=X0, args=ARGS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

pvmismatch-example.ipynb

solar-convex-optimum.png

Note: finding the max power point may be an non-convex optimization problem, therefore to guarantee convergence, find the index of
the max power point by calling numpy.argmax(power). If the spacing of points on the IV curve is too course, then use brentq in a convex trust region around the index of the max power point.