i_sc and v_oc calculations

[Quentin Liebsch]

Hi all,

I am very new to pvlib and currently working on a project where I have to simulate a pv module and right know my results are the dc output (plotted) with values for v_mp, i_mp and p_mp on the y-axis and the time on the x-axis.

I am using pvgis data and calculated the irradiance and cell temperature. Afterwards I use the pvlib.ivtools.sdm.fit_cec_sam function to get the values for: I_L_ref, I_o_ref, R_s, R_sh_ref, a_ref and Adjust and use them to plot my dc output as already written above. 

I was wondering if there is an easy way to also plot the i_sc and v_oc over the time on the x-axis. I also saw the singlediode example for modeling I-V- curves on the pvlib website where they use the desoto function, but they only do that for specific cases with the irradiance and cell temperature. I would like to do that hourly (with the accordingly irradiance and cell temperature) over a specific period to plot the  i_sc and v_oc.

Does anyone have experience with that? (I hope I was being specific enough 😅 )  It would be so helpful! Thank you!

Regards,
Quentin

[cwh...@sandia.gov]

Quentin,

The v_mp, i_mp and p_mp results should be columns in a pandas Dataframe. That same Dataframe should contain i_sc and v_oc.

Can you share your code? 

Cliff

[Quentin Liebsch]

Thank you for the quick response. :) Yes, of course. But I am going to bed soon. I can answer you tomorrow morning again.



import pvlib

from pvlib.location import Location
from pvlib.pvsystem import PVSystem
import pandas as pd
import matplotlib.pyplot as plt


celltype = 'monoSi'
pdc0 = 415
v_mp = 44.1 #Spannung im Maximum Power Point
i_mp = 9.08 #Stromstärke im Maximum Power Point
v_oc = 53.4 #open circuit voltage, wo Linie die x-Achse schneidet I(U)-Graph, für I=0, Leerlaufspannung
i_sc = 9.60 #short circuit current, wo Linie die y-Achse scheidet, für U=0, Kurzschlussstrom
alpha_sc = 0.0005 * i_sc      #Temperaturkoeffizienten für I in %/°C
beta_voc = -0.0029 * v_oc     #Temperaturkoeffizienten für U in %/°C
gamma_pdc = -0.37           #Temperaturkoeffizienten für P in %/°C, muss nicht durch 100 geteilt werden, da der prozentwert abgebildet wird
cells_in_series = 6*27
temp_ref = 25

location = Location(latitude=51.31322846025836, longitude=12.372779830086992,
                    tz='Europe/Berlin', altitude=135, name='NieperBau')

surface_tilt=45
surface_azimuth=0

start="2020-01-01 00:00"
end="2020-01-07 23:00"

poa_data_2020 = pd.read_csv("poa_data_2020_io.csv", index_col=0)
poa_data_2020.index = pd.date_range(start="2020-01-01 00:00",
                                    periods=len(poa_data_2020.index),
                                    freq="h")

poa_data = poa_data_2020[start:end]

solarpos = location.get_solarposition(times=pd.date_range(start=start, end=end
                                                          , freq="h"))
aoi = pvlib.irradiance.aoi(
    surface_tilt, surface_azimuth, solarpos.apparent_zenith, solarpos.azimuth)
iam = pvlib.iam.ashrae(aoi)
effective_irradiance = poa_data["poa_direct"] * iam + poa_data["poa_diffuse"]

temp_cell = pvlib.temperature.faiman(poa_data["poa_global"], poa_data["temp_air"],
                                     poa_data["wind_speed"])



I_L_ref, I_o_ref, R_s, R_sh_ref, a_ref, Adjust = pvlib.ivtools.sdm.fit_cec_sam(
    celltype = celltype,
    v_mp = v_mp,              
    i_mp = i_mp,            
    v_oc = v_oc,             
    i_sc = i_sc,          
    alpha_sc =  alpha_sc,      
    beta_voc = beta_voc,     
    gamma_pmp = gamma_pdc,          
    cells_in_series = cells_in_series,
    temp_ref = temp_ref)


#values for reference conditions from above are now applied to the acutual irradiance data and corresponding cell temperatures
cec_params = pvlib.pvsystem.calcparams_cec(effective_irradiance,  #calculated before
                              temp_cell,                                #calculated from the faimann model
                              alpha_sc,
                              a_ref,
                              I_L_ref,
                              I_o_ref,                          #the rest is calculated by the fitting function above
                              R_sh_ref,
                              R_s,
                              Adjust)


mpp = pvlib.pvsystem.max_power_point(*cec_params,      
                                     method='newton')

# mpp.plot(figsize=(16,9))
# # plt.show()

system = PVSystem(modules_per_string=1, strings_per_inverter=1) #definition des systems, um auch mehrere Zellen kombinieren zu können und auszurechnen
dc_scaled = system.scale_voltage_current_power(mpp) #dc output

dc_scaled.plot(figsize=(16,9))
plt.show()

[cwh...@sandia.gov]

You can get i_sc and v_oc if you replace this line

mpp = pvlib.pvsystem.max_power_point(*cec_params,      
                                     method='newton')

with

dc_results = pvlib.pvsystem.singlediode(*cecparams, method='newton')

mpp = dc_results['p_mp']

i_sc = dc_results['i_sc']

v_oc = dc_results['v_oc']

etc.

Cheers,

Cliff

[Quentin Liebsch]

Thanks for the advice! I will try that out.

Quentin

[Quentin Liebsch]

Hi, 

it's me again and I have to ask for help again. 


now I am getting this error message:

runfile('C:/Users/quent/.spyder-py3/untitled1.py', wdir='C:/Users/quent/.spyder-py3')
Traceback (most recent call last):

  File E:\Anaconda\lib\site-packages\pandas\core\generic.py:554 in _get_axis_number
    return cls._AXIS_TO_AXIS_NUMBER[axis]

KeyError: 1


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  File E:\Anaconda\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\users\quent\.spyder-py3\untitled1.py:91

    dc_scaled = system.scale_voltage_current_power(mpp) #dc output

  File E:\Anaconda\lib\site-packages\pvlib\pvsystem.py:64 in f
    x = func(*args, **kwargs)

  File E:\Anaconda\lib\site-packages\pvlib\pvsystem.py:1063 in scale_voltage_current_power
    return tuple(

  File E:\Anaconda\lib\site-packages\pvlib\pvsystem.py:1064 in <genexpr>
    scale_voltage_current_power(data,

  File E:\Anaconda\lib\site-packages\pvlib\pvsystem.py:3167 in scale_voltage_current_power
    voltage_df = data.filter(voltage_keys, axis=1) * voltage

  File E:\Anaconda\lib\site-packages\pandas\core\generic.py:5450 in filter
    labels = self._get_axis(axis)

  File E:\Anaconda\lib\site-packages\pandas\core\generic.py:566 in _get_axis
    axis_number = self._get_axis_number(axis)

  File E:\Anaconda\lib\site-packages\pandas\core\generic.py:556 in _get_axis_number
    raise ValueError(f"No axis named {axis} for object type {cls.__name__}")

ValueError: No axis named 1 for object type Series


My understanding is that it has something to do with the mpp being a series, but the function system.scale_voltage_current_power is expecting a dataframe. Do I also have to rewrite my last line for it to work?

[cwh...@sandia.gov]

Hi Quentin,

I think you can pass dc_results (the output of pvlib.pvsystem.singlediode) to system.scale_voltage_current_power, the scale function will find the right columns in the DataFrame. It will return a DataFrame with the voltages, currents and power scaled.

Then you can extract i_sc, p_mp, etc. from the returned DataFrame, and the values will be for the system rather than for a single modules.

Cheers,

Cliff

[cwh...@sandia.gov]

I would double check the values of temp_cell and compare with the ambient air temperature from your weather data (call it temp_air). I would expect roughly 30C greater than temp_air in the middle of a day.

I would also plot effective_irradiance vs. dc_scaled['p_mp']. If you don't see a roughly linear pattern, then there's something still wrong in the calculations.

What quantities are in the weather data? The code that gets from the poa_data_2018 being read to effective_irradiance isn't shown.

Is the datetimeindex for the input weather data set to the correct timezone? I don't see it localized in the code posted above.

Cliff

On Thursday, September 7, 2023 at 10:03:16 AM UTC-6 awesome...@gmail.com wrote:

Hi Cliff,

thanks for the advice. Yes, I was able to use the outputs of the singlediode function. I didnt really understand what it did before but now I do. :)

But I am still having some problems with my plots. I don't understand why my diagrams have this unusual peak from mid April to September. I tried this for 3 different years now. Do you have an idea what the cause could be?

my code looks like this now:

**********************************************************************

import pvlib

from pvlib.location import Location
from pvlib.pvsystem import PVSystem
import pandas as pd
import matplotlib.pyplot as plt

import seaborn as sns

celltype = 'monoSi'
pdc0 = 415
v_mp = 44.1 #Spannung im Maximum Power Point
i_mp = 9.08 #Stromstärke im Maximum Power Point
v_oc = 53.4 #open circuit voltage, wo Linie die x-Achse schneidet I(U)-Graph, für I=0, Leerlaufspannung
i_sc = 9.60 #short circuit current, wo Linie die y-Achse scheidet, für U=0, Kurzschlussstrom
alpha_sc = 0.0005 * i_sc      #Temperaturkoeffizienten für I in %/°C
beta_voc = -0.0029 * v_oc     #Temperaturkoeffizienten für U in %/°C
gamma_pdc = -0.37           #Temperaturkoeffizienten für P in %/°C, muss nicht durch 100 geteilt werden, da der prozentwert abgebildet wird
cells_in_series = 6*27
temp_ref = 25

location = Location(latitude=51.31322846025836, longitude=12.372779830086992,
                    tz='Europe/Berlin', altitude=135, name='NieperBau')

surface_tilt=45
surface_azimuth=0

start="2018-01-01 00:00"
end="2018-12-31 23:00"

poa_data_2018 = pd.read_csv("poa_data_2018_io.csv", index_col=0)
poa_data_2018.index = pd.date_range(start="2018-01-01 00:00",
                                    periods=len(poa_data_2018.index),
                                    freq="h")

poa_data = poa_data_2018[start:end]

mpp = pvlib.pvsystem.max_power_point(*cec_params,   #by typing *cec_params you unpack the tuple to use it
                                     method='newton')



# Define empty lists to store results
v_oc_curve_time = []
i_sc_curve_time = []

IL, I0, Rs, Rsh, nNsVth = cec_params


#dann werden die berechneten werte in einer variable (params) gespeichert
params = {                              
    'photocurrent': IL,
    'saturation_current': I0,
    'resistance_series': Rs,
    'resistance_shunt': Rsh,
    'nNsVth': nNsVth
}


curve_info = pvlib.pvsystem.singlediode(method='lambertw', **params)

#Erstellen der v_oc und i_sc Kurven über der Zeit
v_oc_curve_time_df = pd.DataFrame({'v_oc_curve_time': curve_info['v_oc']})
i_sc_curve_time_df = pd.DataFrame({'i_sc_curve_time': curve_info['i_sc']})

# Plot v_oc_curve
v_oc_curve_time_df.plot(figsize=(16, 9), legend=False)
plt.xlabel('Time')
plt.ylabel('Open Circuit Voltage (V)')
plt.title('Open Circuit Voltage Curve')
plt.show()

# Plot i_sc_curve
i_sc_curve_time_df.plot(figsize=(16, 9), legend=False)
plt.xlabel('Time')
plt.ylabel('Short Circuit Current (A)')
plt.title('Short Circuit Current Curve')

plt.show()


system = PVSystem(modules_per_string=1, strings_per_inverter=1) #definition des systems, um auch mehrere Zellen kombinieren zu können und auszurechnen
dc_scaled = system.scale_voltage_current_power(mpp) #dc output

dc_scaled.plot(figsize=(16,9))
plt.show()

************************************************

Thanks in advance!

[kevina...@gmail.com]

Perhaps the surface_azimuth is incorrect?  surface_azimuth=0 means facing north to pvlib.  Such extreme seasonality is not unexpected for a steeply tilted north-facing system at latitude +51, since the module will only receive direct irradiance (aoi < 90) in the summer.

If you intended the system to face south, set surface_azimuth to 180 instead of zero. 

Cheers,

Kevin

[Quentin Liebsch]

Hi,

sorry I deleted my message because I thought it was the wrong code. But it was just shortened, my bad haha.

thank you, I will check that data. 

In the weather data are 8760 values for poa direct, poa sky diffuse, poa ground diffuse, solar elevation, air temperature, wind speed, int, poa diffuse, poa global.


I preprocessed my wether data here:
**************************************************************
import pandas as pd
import pvlib

poa_data_2018, meta, inputs = pvlib.iotools.get_pvgis_hourly(       
    latitude=51.31322846025836, longitude=12.372779830086992,       
    start=2018, end=2018,
    raddatabase="PVGIS-SARAH2", components=True, surface_tilt=45, surface_azimuth=0, 
    outputformat='json', usehorizon=True, userhorizon=None,
    pvcalculation=False, peakpower=None, pvtechchoice='crystSi',
    mountingplace='free', loss=0, trackingtype=0, optimal_surface_tilt=False,
    optimalangles=False, url='https://re.jrc.ec.europa.eu/api/v5_2/',      

   map_variables=True, timeout=30)

poa_data_2018['poa_diffuse'] = poa_data_2018['poa_sky_diffuse'] + poa_data_2018['poa_ground_diffuse']
poa_data_2018['poa_global'] = poa_data_2018['poa_diffuse'] + poa_data_2018['poa_direct']

# poa_data_2020.index = pd.to_datetime(poa_data_2020.index, format="%Y%m%d:%H%M")
print(poa_data_2018)

poa_data_2018.to_csv("poa_data_2018_io.csv")
*********************************************************************

but you are right I didn't set the timezone, but I did it for the location in the main code. But maybe something went wrong here. 

And you are right with the azimuth angle, I wanted to put 180 for south but I got it mixed up with the 0 that I needed for pvgis. Thanks!


Quentin