Improved max. power calculation in pvlib.singlediode._lambertw

revised_lambertw_mpp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed with new method:  0.581497 s\n",
      "Elapsed with old method:  3.603965 s\n",
      "Elapsed with b88 method:  0.525950 s\n",
      "1.0106759873451665e-10\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Created on Tue Jan 24 12:19:19 2023\n",
    "\n",
    "@author: cliff\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "from scipy.special import lambertw\n",
    "from scipy.optimize import brentq, newton\n",
    "from pvlib.singlediode import (\n",
    "    _lambertw_i_from_v, _lambertw_v_from_i, _lambertw, bishop88_mpp, bishop88_i_from_v, bishop88_v_from_i,\n",
    "    estimate_voc, _prepare_newton_inputs, bishop88)\n",
    "\n",
    "from timeit import timeit\n",
    "\n",
    "\n",
    "# code for secant method to get MP values\n",
    "\n",
    "def w_psi(i, il, io, rs, rsh, a):\n",
    "    gsh = 1. / rsh\n",
    "    with np.errstate(over='ignore'):\n",
    "        argW = (io / (gsh * a) *\n",
    "                np.exp((-i + il + io) /\n",
    "                       (gsh * a)))\n",
    "    lambertwterm = np.array(lambertw(argW).real)\n",
    "\n",
    "    idx_inf = np.logical_not(np.isfinite(lambertwterm))\n",
    "    if np.any(idx_inf):\n",
    "        # Calculate using log(argW) in case argW is really big\n",
    "        logargW = (np.log(io) - np.log(gsh) -\n",
    "                   np.log(a) +\n",
    "                   (-i + il + io) /\n",
    "                   (gsh * a))[idx_inf]\n",
    "\n",
    "        # Three iterations of Newton-Raphson method to solve\n",
    "        #  w+log(w)=logargW. The initial guess is w=logargW. Where direct\n",
    "        #  evaluation (above) results in NaN from overflow, 3 iterations\n",
    "        #  of Newton's method gives approximately 8 digits of precision.\n",
    "        w = logargW\n",
    "        for _ in range(0, 6):\n",
    "            w = w * (1. - np.log(w) + logargW) / (1. + w)\n",
    "        lambertwterm[idx_inf] = w\n",
    "    return lambertwterm\n",
    "\n",
    "\n",
    "def f(i, il, io, rs, rsh, a):\n",
    "    wma = w_psi(i, il, io, rs, rsh, a)\n",
    "    return (il + io - i) * rsh - i * rs - a * wma\n",
    "\n",
    "\n",
    "def fprime(i, il, io, rs, rsh, a):\n",
    "    wma = w_psi(i, il, io, rs, rsh, a)\n",
    "    return -rs - rsh / (1 + wma)\n",
    "\n",
    "\n",
    "def fprime2(i, il, io, rs, rsh, a):\n",
    "    wma = w_psi(i, il, io, rs, rsh, a)\n",
    "    return -rsh**2. / a * wma / (1 + wma)**3.\n",
    "\n",
    "\n",
    "def imp_est(i, il, io, rs, rsh, a):\n",
    "    wma = w_psi(i, il, io, rs, rsh, a)\n",
    "    f = (il + io - i) * rsh - i * rs - a * wma\n",
    "    fprime = -rs - rsh / (1 + wma)\n",
    "    return -f / fprime\n",
    "\n",
    "\n",
    "def imp_zero(i, il, io, rs, rsh, a):\n",
    "    return imp_est(i, il, io, rs, rsh, a) - i\n",
    "\n",
    "\n",
    "def imp_zero_prime(i, il, io, rs, rsh, a):\n",
    "    wma = w_psi(i, il, io, rs, rsh, a)\n",
    "    f = (il + io - i) * rsh - i * rs - a * wma\n",
    "    fprime = -rs - rsh / (1 + wma)\n",
    "    fprime2 = -rsh**2. / a * wma / (1 + wma)**3.\n",
    "    return f / fprime**2. * fprime2 - 2.\n",
    "\n",
    "\n",
    "# alternative to pvlib.singlediode._lambertw, replaces GSS with newton + analytic method\n",
    "\n",
    "def _lambertw_new(photocurrent, saturation_current, resistance_series,\n",
    "              resistance_shunt, nNsVth, ivcurve_pnts=None):\n",
    "    # Compute short circuit current\n",
    "    i_sc = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth, 0.,\n",
    "                              saturation_current, photocurrent)\n",
    "\n",
    "    # Compute open circuit voltage\n",
    "    v_oc = _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth, 0.,\n",
    "                              saturation_current, photocurrent)\n",
    "\n",
    "    params = (photocurrent, saturation_current, resistance_series,\n",
    "              resistance_shunt, nNsVth)\n",
    "    \n",
    "    # Compute maximum power point quantities\n",
    "    imp_est = 0.8 * photocurrent\n",
    "    args, i0 = _prepare_newton_inputs((), params, imp_est)\n",
    "    i_mp = newton(func=imp_zero, x0=i0, fprime=imp_zero_prime,\n",
    "                 args=args)\n",
    "    v_mp = _lambertw_v_from_i(resistance_shunt, resistance_series, nNsVth, i_mp,\n",
    "                              saturation_current, photocurrent)\n",
    "    p_mp = i_mp * v_mp\n",
    "\n",
    "    # Find Ix and Ixx using Lambert W\n",
    "    i_x = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth,\n",
    "                             0.5 * v_oc, saturation_current, photocurrent)\n",
    "\n",
    "    i_xx = _lambertw_i_from_v(resistance_shunt, resistance_series, nNsVth,\n",
    "                              0.5 * (v_oc + v_mp), saturation_current,\n",
    "                              photocurrent)\n",
    "\n",
    "    out = (i_sc, v_oc, i_mp, v_mp, p_mp, i_x, i_xx)\n",
    "\n",
    "    # create ivcurve\n",
    "    if ivcurve_pnts:\n",
    "        ivcurve_v = (np.asarray(v_oc)[..., np.newaxis] *\n",
    "                     np.linspace(0, 1, ivcurve_pnts))\n",
    "\n",
    "        ivcurve_i = _lambertw_i_from_v(resistance_shunt, resistance_series,\n",
    "                                       nNsVth, ivcurve_v.T, saturation_current,\n",
    "                                       photocurrent).T\n",
    "\n",
    "        out += (ivcurve_i, ivcurve_v)\n",
    "\n",
    "    return out\n",
    "\n",
    "\n",
    "def new_way(params):\n",
    "    return _lambertw_new(*params)\n",
    "\n",
    "    \n",
    "def old_way(params):\n",
    "    return _lambertw(*params)\n",
    "\n",
    "\n",
    "def b88(params):\n",
    "    # for timing, need to also compute isc, voc, i_x and i_xx\n",
    "    isc = bishop88_i_from_v(0., *params)\n",
    "    voc = bishop88_v_from_i(0., *params)\n",
    "    imp, vmp, pmp = bishop88_mpp(*params)\n",
    "    ix = bishop88_i_from_v(vmp/2., *params)\n",
    "    ixx = bishop88_i_from_v((voc + vmp)/2., *params)\n",
    "    return imp, vmp, pmp\n",
    "\n",
    "\n",
    "base_params = (1., 5.e-9, 0.5, 2000., 72 * 1.1 * 0.025)\n",
    "nsamples = 10000\n",
    "\n",
    "il = 9 * np.random.rand(nsamples) + 1.  # 1.- 10. A\n",
    "io = 10**(-9 + 3. * np.random.rand(nsamples))  # 1e-9 to 1e-6 A\n",
    "rs = 5 * np.random.rand(nsamples) + 0.05  # 0.05 to 5.05 Ohm\n",
    "rsh = 10**(2 + 2 * np.random.rand(nsamples))  # 100 to 10000 Ohm\n",
    "n = 1 + 0.7 * np.random.rand(nsamples)  #  1.0 to 1.7\n",
    "nNsVth = 72 * n * 0.025\n",
    "params = (il, io, rs, rsh, nNsVth)\n",
    "\n",
    "elapsed_new = timeit('new_way(params)', number=10, globals=globals())\n",
    "elapsed_old = timeit('old_way(params)', number=10, globals=globals())\n",
    "elapsed_b88 = timeit('b88(params)', number=10, globals=globals())\n",
    "\n",
    "print('Elapsed with new method: %9.6f s' % elapsed_new)\n",
    "print('Elapsed with old method: %9.6f s' % elapsed_old)\n",
    "print('Elapsed with b88 method: %9.6f s' % elapsed_b88)\n",
    "\n",
    "# check for agreement\n",
    "b88_result = bishop88_mpp(*params)\n",
    "result = _lambertw_new(*params)\n",
    "print(np.abs(result[4] - b88_result[2]).max())\n"
   ]
  },
  {
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}