Convergence issues with gen_coeffs

[mikofski]

Several folks have had convergence issues with gen_coeffs, but there are some tricks that can help. Under the hood it just uses scipy.optimize.root, which depends strongly on the initial guess. It also has interesting criteria for exiting the solver. Therefore it's important to look at the solver message and the residuals to determine why scipy.optimize.root decided to quit. Sometimes it's not making progress on the solution, but it actually has found a good fit in the least squares sense. The message in this case is usually:

The iteration is not making good progress, as measured by the improvement from the last ten iterations.

The solver solution shows the residuals in sol.fun. So if the solver quits due to lack of good progress, and the values of sol.fun are all extremely small, that's probably the best solution you can get.

Another trick I've found that seems to work sometimes is to restart the solver with the last guess. The 2-diode equations are "stiff" so often one of the guess "blows up" and needs to be reset to a reasonable value before continuing.

Take for example this Canadian Solar CS6X-300M module:

Name	Vintage	Area	Material	Cells in Series	Parallel Strings	Isco	Voco	Impo	Vmpo	Aisc
Canadian Solar CS6X-300M [2013]	2013	1.91	c-Si	72	1	8.6388	43.5918	8.1359	34.9531	0.0005

import numpy as np
import pvlib
from pvmismatch import *
from pvmismatch.contrib import gen_coeffs

# get the sandia modules from SAM using pvlib
sandia_mod = pvlib.pvsystem.retrieve_sam('SandiaMod')
# search for the Canadian Solar module
sandia_mod.columns[sandia_mod.T.index.str.startswith('Canadian')]
# this returns 2 modules:
# Index(['Canadian_Solar_CS5P_220M___2009_', 'Canadian_Solar_CS6X_300M__2013_'], dtype='object')

pvmod = sandia_mod.columns[sandia_mod.T.index.str.startswith('Canadian')][1]
pvmod = sandia_mod[pvmod]

args = (pvmod.Isco, pvmod.Voco, pvmod.Impo, pvmod.Vmpo, pvmod.Cells_in_Series, pvmod.Parallel_Strings, 25.0)
x, sol = gen_coeffs.gen_two_diode(*args)
sol  # solution results
#    fjac: array([[-7.66371060e-01, -4.12098681e-02, -6.41074992e-01,
#        -1.11385313e-09],
#       [ 6.37989585e-01, -1.65538611e-01, -7.52041393e-01,
#        -2.63306547e-05],
#       [ 7.51311392e-02,  9.85341900e-01, -1.53155468e-01,
#         2.32875163e-04],
#       [ 6.98346451e-07,  2.33820448e-04, -1.58636487e-05,
#        -9.99999973e-01]])
#     fun: array([ 1.97339520e-08,  6.92527653e-08,  7.86377008e-09, -2.31607963e-06])
# message: 'The iteration is not making good progress, as measured by the \n  improvement from the last ten iterations.'
#    nfev: 223
#    njev: 9
#     qtf: array([-2.28610040e-08, -5.32559619e-09,  6.30173323e-08,  2.31701635e-06])
#       r: array([ 1.10432017e+01,  2.61218295e-01,  2.05879162e+02, -4.37501467e-03,
#        1.04169191e-01,  2.44093888e+02, -1.93746667e-05,  2.69645844e+01,
#        2.97703914e-03,  7.57173785e-07])
#  status: 5
# success: False
#       x: array([-21.42656746, -13.75511821,   0.07434147,   6.11722022])

x  # b/c the solver "failed" this is the initial guess x0 which defaults to pvcell.ISAT1_T0, pvcell.ISAT2_T0, pvcell.RS, pvcell.RSH
# (2.28618816125344e-11,
#  1.117455042372326e-06,
#  0.004267236774264931,
#  10.01226369025448)

So to restart the solution we need the last guess, which was transformed to log-space and parabolic-space to help it solve better, b/c we expect Isat1 and Isat2 to be positive and to change by an order of magnitude or more and Rs and Rsh to be positive.

def last_guess(sol):
    isat1 = np.exp(sol.x[0])
    isat2 = np.exp(sol.x[1])
    rs = sol.x[2] ** 2.0
    rsh = sol.x[3] ** 2.0
    return isat1, isat2, rs, rsh

x = last_guess(sol)
# (4.949484718740044e-10,
#  1.0622532408228556e-06,
#  0.0055266539292632725,
#  37.420383187251204)
# these don't actually seem that bad, so no changes required

# restart the solver
x, sol = gen_coeffs.gen_two_diode(*args, x0=x)
#    fjac: array([[-7.60870868e-01, -4.58108789e-02, -6.47284239e-01,
#         3.72235582e-07],
#       [ 6.40027744e-01, -2.17398354e-01, -7.36954844e-01,
#        -2.37286186e-05],
#       [ 1.06957980e-01,  9.75007332e-01, -1.94732309e-01,
#         1.43209394e-04],
#       [-1.52809906e-07,  1.44771721e-04, -1.06415178e-05,
#        -9.99999989e-01]])
#     fun: array([ 1.97339520e-08,  6.92527653e-08,  7.86377008e-09, -2.31607963e-06])
# message: 'The iteration is not making good progress, as measured by the \n  improvement from the last ten iterations.'
#    nfev: 11
#    njev: 2
#     qtf: array([-2.32784758e-08, -8.16544653e-09,  6.77696431e-08,  2.31608954e-06])
#       r: array([ 1.13359731e+01,  2.76160043e-01,  2.15419870e+02, -6.66722799e-03,
#        1.04100761e-01,  2.43542797e+02, -2.39314749e-04,  3.97725352e+01,
#        4.44961882e-03,  7.68389446e-07])
#  status: 5
# success: False
#       x: array([-21.42656746, -13.75511821,   0.07434147,   6.11722022])

So it's stll not making progress, and in fact there's almost no change from the last run other than the Jacobian is slightly different. we could run it again just in case

x = last_guess(sol)
# (4.949484718740044e-10,
#  1.0622532408228556e-06,
#  0.0055266539292632725,
#  37.420383187251204)
# exactly the same

x, sol = gen_coeffs.gen_two_diode(*args, x0=x)
# this time sol is exactly the same, so it's not going to change unless we change x0

# instead take a look at sol.fun
# array([ 1.97339520e-08,  6.92527653e-08,  7.86377008e-09, -2.31607963e-06])
# these are very small numbers, so let's just try the last guess in a PVcell and see how well it fits?

pvc = pvcell.PVcell(Isat1_T0=x[0], Isat2_T0=x[1], Rs=x[2], Rsh=x[3], Isc0_T0=pvmod.Isco, alpha_Isc=pvmod.Aisc)

pvc.Isc
# 8.6388 [A] ha-ha! this is always exactly the same because it's an argument to PVcell
pvmod.Isco
# 8.6388 [A]

pvc.Voc*72
# 43.59529549470483 [V]
pvmod.Voco
# 43.5918 [V] so off by about 0.01% off, that is really close! 

mpp = np.argmax(pvc.Pcell)  # find the index of the max power point

pvc.Icell[mpp]
# array([8.12227748])
pvmod.Impo
# 8.1359 [A] so this is off by about 0.17%, not zero but not too shabby

pvc.Vcell[mpp]*72
# array([35.01088617])
pvmod.Vmpo
# 34.9531 [V] also about 0.17% off, which is pretty good

# check the fill factor:
(pvmod.Impo * pvmod.Vmpo) / pvmod.Isco / pvmod.Voco
# 0.7551497439448178
(pvc.Icell[mpp] * pvc.Vcell[mpp]) / pvc.Isc / pvc.Voc
# array([0.75507116])

Okay, we found a solution that's pretty good actually, but just for fun, let's still see if we can do better? This is a totally random swag, but b/c the fill factor of the fit is slightly low, let's increase Rsh just to see what happens?

x, sol = gen_coeffs.gen_two_diode(*args, x0=(x[0], x[1], x[2], 100))
#    fjac: array([[-7.60738779e-01, -4.16279681e-02, -6.47721871e-01,
#        -1.77631852e-08],
#       [ 6.39635461e-01, -2.17509958e-01, -7.37262432e-01,
#        -2.46843217e-05],
#       [ 1.10195221e-01,  9.75169989e-01, -1.92094985e-01,
#         1.48542234e-04],
#       [ 5.93189989e-07,  1.50223756e-04, -1.03238897e-05,
#        -9.99999989e-01]])
#     fun: array([ 4.27721605e-10, -7.93960453e-10,  7.63408892e-10, -1.39545171e-08])
# message: 'The number of calls to function has reached maxfev = 500.'
#    nfev: 500
#    njev: 6
#     qtf: array([-7.87560405e-10, -1.13267653e-10, -8.79397959e-10,  1.40880431e-08])
#       r: array([ 1.12788187e+01,  7.63051802e-02,  2.12444826e+02,  1.69125666e-02,
#        2.91168216e-02,  2.45221869e+02,  5.14488430e-04,  4.07251052e+01,
#       -1.15779878e-02, -2.35233507e-06])
#  status: 2
# success: False
#       x: array([-21.41642896, -15.02715388,   0.07464628,  -4.41870535])

This time the solver iterated 500 times, the default max, but was still making good progress so start again!

x0 = last_guess(sol)
# (4.999920303677904e-10,
# 2.977076485134717e-07,
# 0.005572067153316041,
# 19.524956939393345)

x0, sol = gen_coeffs.gen_two_diode(*args, x0)
#    fjac: array([[-7.59355687e-01, -4.01508771e-02, -6.49435792e-01,
#        -3.62165935e-08],
#       [ 6.40821098e-01, -2.19198593e-01, -7.35731131e-01,
#        -2.44857662e-05],
#       [ 1.12815162e-01,  9.74853765e-01, -1.92179228e-01,
#         1.45315821e-04],
#       [ 7.30333726e-07,  1.47030377e-04, -9.88822165e-06,
#        -9.99999989e-01]])
#     fun: array([ 6.01669096e-11, -2.54640753e-11,  1.28311584e-10, -1.15756661e-10])
# message: 'The solution converged.'
#    nfev: 275
#    njev: 2
#     qtf: array([-1.27995811e-10, -5.02620775e-11, -4.27116523e-11,  1.15751691e-10])
#       r: array([ 1.12947984e+01,  7.54656350e-02,  2.13349054e+02, -1.71717846e-02,
#        2.87152855e-02,  2.45323088e+02, -5.18920402e-04,  4.13221785e+01,
#        1.17930254e-02,  2.36371524e-06])
#  status: 1
# success: True
#       x: array([-21.41637331, -15.04143894,   0.07464795,   4.41283666])
# TADA! we solved it! Notice sol.fun is even smaller than it was for the last solution!

# now check the cell characteristics again
pvc0 = pvcell.PVcell(Isat1_T0=x0[0], Isat2_T0=x0[1], Rs=x0[2], Rsh=x0[3], Isc0_T0=pvmod.Isco, alpha_Isc=pvmod.Aisc)

pvc0.Isc
# 8.6388 [A] ha-ha! this is always exactly what you set it to in PVcell()

pvc0.Voc*72
# 43.598482683059295 [V] this really isn't too different from the previous solution
pvc.Voc*72
# 43.59529549470483 [V]
pvmod.Voco
# 43.5918 [V]

mpp0 = np.argmax(pvc0.Pcell)
pvc0.Icell[mpp0]
# array([8.12854724]) this is also very close to the last solution, but the error is 0.09%, about half!
pvc.Icell[mpp]
# array([8.12227748])
In [83]: pvmod.Impo
Out[83]: 8.1359

pvc0.Vcell[mpp0]*72
# array([34.98446279]) again super close to last solution, but the error has been halved to 0.09%!
pvc.Vcell[mpp]*72
# array([35.01088617])
pvmod.Vmpo
# 34.9531 [V]

(pvc0.Icell[mpp0] * pvc0.Vcell[mpp0]) / pvc0.Isc / pvc0.Voc
# array([0.75502851]) the fill factor actually stayed about the same

Anyway there you have it; a few tricks to improve convergence:

1. check the residuals, sol.fun and solver message, sol.message
2. restart with the last guess, tweaking values if needed
3. use sol.x and un-transform it to get the last guess

   def last_guess(sol):
       isat1 = np.exp(sol.x[0])
       isat2 = np.exp(sol.x[1])
       rs = sol.x[2] ** 2.0
       rsh = sol.x[3] ** 2.0
       return isat1, isat2, rs, rsh

4. create a test PVcell and check Voc, Icell[np.argmax(Pcell)], Vcell[np.argmax(Pcell)], and the fill factor

[mikofski]

here's another example:

"""
Making 2-diode modules from CEC in PVMismatch
"""

from matplotlib import pyplot as plt
import numpy as np
import pvlib
from pvmismatch import *
from pvmismatch.contrib import gen_coeffs

cecmod = pvlib.pvsystem.retrieve_sam('CECMod')
csmods = cecmod.columns[cecmod.T.index.str.startswith('Canadian')]
len(csmods)  # 409 modules in CEC library!
cs6x_300m = cecmod[csmods[264]]  # CS6X-300M Canadian Solar 300W mono-Si module

args = (
    cs6x_300m.I_sc_ref,
    cs6x_300m.V_oc_ref,
    cs6x_300m.I_mp_ref,
    cs6x_300m.V_mp_ref,
    cs6x_300m.N_s,  # number of series cells
    1,              # number of parallel sub-strings
    25.0)           # cell temperature

# try to solve using default coeffs
x, sol = gen_coeffs.gen_two_diode(*args)

# solver fails, so get the last guess before it quit
def last_guess(sol):
    isat1 = np.exp(sol.x[0])
    isat2 = np.exp(sol.x[1])
    rs = sol.x[2] ** 2.0
    rsh = sol.x[3] ** 2.0
    return isat1, isat2, rs, rsh
x = last_guess(sol)

# the series and shunt resistance solver guess are way off, so reset them
# with something reasonable for a 2-diode model
x, sol = gen_coeffs.gen_two_diode(*args, x0=(x[0], x[1], 0.005, 10.0))

# Hooray, it worked! Note that the 1-diode and 2-diode parametres are so
# different! Anyway, let's make a cell and a module to check the solution.
pvc = pvcell.PVcell(
    Isat1_T0=x[0],
    Isat2_T0=x[1],
    Rs=x[2],
    Rsh=x[3],
    Isc0_T0=cs6x_300m.I_sc_ref,
    alpha_Isc=cs6x_300m.alpha_sc)
np.isclose(pvc.Isc, cs6x_300m.I_sc_ref)  # ha-ha, this exact b/c we used it
# open circuit voltage within 1E-3: (45.01267251639085, 45.0)
np.isclose(pvc.Voc*cs6x_300m.N_s, cs6x_300m.V_oc_ref, rtol=1e-3, atol=1e-3)

# get index max power point
mpp = np.argmax(pvc.Pcell)

# max power voltage within 1E-3: (36.50580418834946, 36.5)
np.isclose(pvc.Vcell[mpp][0]*cs6x_300m.N_s, cs6x_300m.V_mp_ref, rtol=1e-3, atol=1e-3)
# max power current within 1E-3: (8.218687568902466, 8.22)
np.isclose(pvc.Icell[mpp][0], cs6x_300m.I_mp_ref, rtol=1e-3, atol=1e-3)

# use pvlib to get the full IV curve using CEC model
params1stc = pvlib.pvsystem.calcparams_cec(effective_irradiance=1000,
    temp_cell=25.0, alpha_sc=cs6x_300m.alpha_sc, a_ref=cs6x_300m.a_ref,
    I_L_ref=cs6x_300m.I_L_ref, I_o_ref=cs6x_300m.I_o_ref, R_sh_ref=cs6x_300m.R_sh_ref,
    R_s=cs6x_300m.R_s, Adjust=cs6x_300m.Adjust)
iv_params1stc = pvlib.pvsystem.singlediode(*params1stc, ivcurve_pnts=100, method='newton')

# use pvmm to get full IV curve using 2-diode model parameters
pvm = pvmodule.PVmodule(cell_pos=pvmodule.STD72, pvcells=[pvc]*72)

# make some comparison plots
pvm.plotMod()  # plot the pvmm module
plt.tight_layout()

# get axes for IV curve
f, ax = plt.gcf(), plt.gca()
ax0 = f.axes[0]
ax0.plot(iv_params1stc['v'], iv_params1stc['i'], '--')
ax0.plot(0, iv_params1stc['i_sc'], '|k')
ax0.plot(iv_params1stc['v_oc'], 0, '|k')
ax0.plot(iv_params1stc['v_mp'], iv_params1stc['i_mp'], '|k')
ax0.set_ylim([0, 10])

ax0.plot(0, pvm.Isc.mean(), '_k')
ax0.plot(pvm.Voc.sum(), 0, '|k')

mpp = np.argmax(pvm.Pmod)
ax0.plot(pvm.Vcell[mpp], mpp.Icell[mpp], '_k')
ax0.plot(pvm.Vmod[mpp], pvm.Imod[mpp], '_k')

iv_params1stc['p'] = iv_params1stc['v'] * iv_params1stc['i']
ax1.plot(iv_params1stc['v'], iv_params1stc['p'], '--')

ax1.plot(iv_params1stc['v_mp'], iv_params1stc['p_mp'], '|k')
ax1.plot(pvm.Vmod[mpp], pvm.Pmod[mpp], '_k')

[beingwei]

I used the method of this code “def last_guess(sol):
isat1 = np.exp(sol.x[0])
isat2 = np.exp(sol.x[1])
rs = sol.x[2] ** 2.0
rsh = sol.x[3] ** 2.0
return isat1, isat2, rs, rsh
x = last_guess(sol)”, resulting in the following error：
D:\Program Files\Python39\lib\site-packages\pvmismatch\pvmismatch_lib\pvconstants.py:58: RuntimeWarning: divide by zero encountered in true_divide
yright = fp[-1] + (xright - xp[-1]) / (xp[-2] - xp[-1]) * (fp[-2] - fp[-1])
D:\Program Files\Python39\lib\site-packages\pvmismatch\pvmismatch_lib\pvconstants.py:167: RuntimeWarning: invalid value encountered in multiply
Vforward = Vmax * self.pts
D:\Program Files\Python39\lib\site-packages\pvmismatch\pvmismatch_lib\pvconstants.py:53: RuntimeWarning: invalid value encountered in double_scalars
yleft = fp[0] + (xleft - xp[0]) / (xp[1] - xp[0]) * (fp[1] - fp[0])

The final error is as follows：
pvsys = pvm.PVsystem(numberStrs=1, numberMods=1, pvmods=pvmod)
File "D:\Program Files\Python39\lib\site-packages\pvmismatch\pvmismatch_lib\pvsystem.py", line 62, in init
self.update()
File "D:\Program Files\Python39\lib\site-packages\pvmismatch\pvmismatch_lib\pvsystem.py", line 68, in update
self.Isc, self.Voc, self.FF, self.eff) = self.calcMPP_IscVocFFeff()
File "D:\Program Files\Python39\lib\site-packages\pvmismatch\pvmismatch_lib\pvsystem.py", line 108, in calcMPP_IscVocFFeff
Vmp = (-Pv[0] * np.diff(Vmid, axis=0) / np.diff(Pv, axis=0) + Vmid[0]).item()
IndexError: index 0 is out of bounds for axis 0 with size 0

I didn't find the reason

[mikofski]

It’s really hard to know what the problem is without the exact code you are trying to run. Maybe create a gist with your code and paste the link here? Make sure that your code has whatever coefficients you’re using. Also what are you trying to achieve? Do you have a module data sheet that you want to generate 2-diode coefficients for?