From 0d47ed030ec504249981816983fe02b390558ab7 Mon Sep 17 00:00:00 2001 From: Darice Date: Sun, 9 Feb 2025 06:22:16 -0700 Subject: [PATCH 01/17] create file --- files/SixParSolve.py | 89 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 89 insertions(+) create mode 100644 files/SixParSolve.py diff --git a/files/SixParSolve.py b/files/SixParSolve.py new file mode 100644 index 0000000..257824b --- /dev/null +++ b/files/SixParSolve.py @@ -0,0 +1,89 @@ +import pyomo.environ as pyo +import pandas as pd +import numpy as np + +m = pyo.ConcreteModel() + +m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.bVoc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.aIsc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.gPmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + +m.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15)) +m.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20)) +m.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7)) # might need to scale +m.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75)) +m.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6)) +m.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100)) + +m.Tc = pyo.RangeSet(-10 + 273.15, 50 + 273.15, 3) +m.nTc = pyo.Param(len(m.Tc)) +m.V_Tc = pyo.Var(m.Tc, domain=pyo.NonNegativeReals) +m.I_Tc = pyo.Var(m.Tc, domain=pyo.NonNegativeReals) + +m.f0 = pyo.Constraint(expr=m.Il - m.Io*( pyo.exp( m.Isc*m.Rs / m.a ) - 1 ) - m.Isc*m.Rs/m.Rsh - m.Isc == 0) +m.f1 = pyo.Constraint(expr=m.Io*( pyo.exp( m.Voc/m.a ) - 1 ) + m.Voc/m.Rsh -m.Il == 0) +m.f2 = pyo.Constraint(expr=m.Il - m.Io*( pyo.exp( (m.Vmp + m.Imp*m.Rs) / m.a ) - 1 ) - (m.Vmp + m.Imp*m.Rs)/m.Rsh - m.Imp == 0) +m.f3 = pyo.Constraint(expr=m.Imp - m.Vmp*( + ( m.Io/m.a*pyo.exp( (m.Vmp + m.Imp*m.Rs)/m.a ) + 1/m.Rsh ) + /( 1 + m.Io*m.Rs/m.a*pyo.exp( (m.Vmp + m.Imp*m.Rs)/m.a ) + m.Rs/m.Rsh ) ) == 0) + +m.dT = pyo.Param(5) + +m.aT = pyo.Expression(expr=m.a*(m.Tref+m.dT)/m.Tref) +m.VocT = pyo.Expression(expr=m.bVoc*(1+m.Adj/100.0)*m.dT + m.Voc) +m.Eg = pyo.Expression(expr=(1-0.0002677*m.dT)*m.Egref) +m.IoT = pyo.Expression(expr=m.Io*( (m.Tref+m.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.Eg/(m.Tref+m.dT)))) + +m.f4 = pyo.Expression(expr=m.Il+m.aIsc*(1-m.Adj/100)*m.dT - m.IoT*(pyo.exp( m.VocT/m.aT ) - 1 ) - m.VocT/m.Rsh) + +m.a_Tc = pyo.Expression(m.Tc, rule=lambda t: m.a * m.Tc[t] / m.Tref) +m.Io_Tc = pyo.Expression(m.Tc, rule=lambda t: m.Io* ( m.Tc[t]/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.Eg/m.Tc[t]))) +m.Il_Tc = pyo.Expression(m.Tc, rule=lambda t: m.Il + (m.aIsc*(1-m.Adj/100))*(m.Tc[t]-m.Tref)) +m.f_5 = pyo.Constraint(m.Tc, rule=lambda t: m.V_Tc[t] *( m.Io_Tc[t]/m.a_Tc[t]*pyo.exp( (m.V_Tc[t]+m.I_Tc[t]*m.Rs)/m.a_Tc[t] ) + 1/m.Rsh ) + / ( 1 + m.Rs/m.Rsh + m.Io_Tc[t]*m.Rs/m.a_Tc[t]*pyo.exp( (m.V_Tc[t]+m.I_Tc[t]*m.Rs)/m.a_Tc[t] ) ) - m.I_Tc[t] == 0) +m.f_6 = pyo.Constraint(m.Tc, rule=lambda t: m.Il_Tc[t] - m.Io_Tc[t]*(pyo.exp( (m.V_Tc[t]+m.I_Tc[t]*m.Rs)/m.a_Tc[t] ) - 1) - (m.V_Tc[t] + m.I_Tc[t]*m.Rs)/m.Rsh - m.I_Tc[t]) +m.Pmax_Tc = pyo.Expression(m.Tc, rule=lambda t: m.V_Tc[t] * m.I_Tc[t]) + +def gamma_rule(m, t): + if t >= m.nTc: + return pyo.Constraint.Skip + return (m.Pmax_Tc[t]-m.Pmax_Tc[t-1])*100/(m.Vmp*m.Imp*(m.Tc[t]-m.Tc[t-1])) + +m.gamma_Tc = pyo.Expression(m.Tc, rule=gamma_rule) +m.gamma = pyo.Expression(expr=sum(m.gamma_Tc[t] for t in m.Tc) / m.nTc) + +m.f_7 = pyo.Constraint(expr=m.gamma - m.gPmp == 0) + +# Sanity checks on current + +m.Imp = pyo.Var(domain=pyo.NonNegativeReals) +m.f_8 = pyo.Constraint(expr=m.Il - m.Imp - m.Io * (pyo.exp((m.Vmp + m.Imp * m.Rs) / m.a) - 1.0) - (m.Vmp + m.Imp * m.Rs) / m.Rsh) + +m.Ioc = pyo.Var(domain=pyo.NonNegativeReals) +m.f_9 = pyo.Constraint(expr=m.Il - m.Ioc - m.Io * (pyo.exp((m.Voc + m.Ioc * m.Rs) / m.a) - 1.0) - (m.Voc + m.Ioc * m.Rs) / m.Rsh) + +# Sanity checks on slope +v_range = np.linspace(0.015 * m.Voc, 0.98 * m.Voc, 100) + +filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" +with open(filename, 'r') as f: + +if Imp < Isc: + raise ValueError + +Vmp.set_value() +Imp.set_value() +Voc.set_value() +Isc.set_value() +bVoc.set_value() +aIsc.set_value() +gPmp.set_value() +Egref.set_value() +Tref.set_value() + From b389712762ec6d257d65060b3e050a711ab48127 Mon Sep 17 00:00:00 2001 From: Darice Date: Mon, 10 Feb 2025 14:32:24 -0700 Subject: [PATCH 02/17] split into blocks --- files/SixParSolve.py | 177 ++++++++++++++++++++++++++----------------- 1 file changed, 107 insertions(+), 70 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 257824b..0dcff7d 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -1,89 +1,126 @@ import pyomo.environ as pyo import pandas as pd import numpy as np +import idaes.logger as idaeslog +from pyomo.util.infeasible import log_infeasible_constraints, log_infeasible_bounds, log_close_to_bounds -m = pyo.ConcreteModel() +solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties") + +model = pyo.ConcreteModel() +model.solver = pyo.Block() +m = model.solver m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.bVoc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.aIsc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.gPmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) +m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) +m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) +m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) +m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15)) -m.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20)) -m.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7)) # might need to scale -m.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75)) -m.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6)) -m.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100)) - -m.Tc = pyo.RangeSet(-10 + 273.15, 50 + 273.15, 3) -m.nTc = pyo.Param(len(m.Tc)) -m.V_Tc = pyo.Var(m.Tc, domain=pyo.NonNegativeReals) -m.I_Tc = pyo.Var(m.Tc, domain=pyo.NonNegativeReals) - -m.f0 = pyo.Constraint(expr=m.Il - m.Io*( pyo.exp( m.Isc*m.Rs / m.a ) - 1 ) - m.Isc*m.Rs/m.Rsh - m.Isc == 0) -m.f1 = pyo.Constraint(expr=m.Io*( pyo.exp( m.Voc/m.a ) - 1 ) + m.Voc/m.Rsh -m.Il == 0) -m.f2 = pyo.Constraint(expr=m.Il - m.Io*( pyo.exp( (m.Vmp + m.Imp*m.Rs) / m.a ) - 1 ) - (m.Vmp + m.Imp*m.Rs)/m.Rsh - m.Imp == 0) -m.f3 = pyo.Constraint(expr=m.Imp - m.Vmp*( - ( m.Io/m.a*pyo.exp( (m.Vmp + m.Imp*m.Rs)/m.a ) + 1/m.Rsh ) - /( 1 + m.Io*m.Rs/m.a*pyo.exp( (m.Vmp + m.Imp*m.Rs)/m.a ) + m.Rs/m.Rsh ) ) == 0) - -m.dT = pyo.Param(5) - -m.aT = pyo.Expression(expr=m.a*(m.Tref+m.dT)/m.Tref) -m.VocT = pyo.Expression(expr=m.bVoc*(1+m.Adj/100.0)*m.dT + m.Voc) -m.Eg = pyo.Expression(expr=(1-0.0002677*m.dT)*m.Egref) -m.IoT = pyo.Expression(expr=m.Io*( (m.Tref+m.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.Eg/(m.Tref+m.dT)))) - -m.f4 = pyo.Expression(expr=m.Il+m.aIsc*(1-m.Adj/100)*m.dT - m.IoT*(pyo.exp( m.VocT/m.aT ) - 1 ) - m.VocT/m.Rsh) - -m.a_Tc = pyo.Expression(m.Tc, rule=lambda t: m.a * m.Tc[t] / m.Tref) -m.Io_Tc = pyo.Expression(m.Tc, rule=lambda t: m.Io* ( m.Tc[t]/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.Eg/m.Tc[t]))) -m.Il_Tc = pyo.Expression(m.Tc, rule=lambda t: m.Il + (m.aIsc*(1-m.Adj/100))*(m.Tc[t]-m.Tref)) -m.f_5 = pyo.Constraint(m.Tc, rule=lambda t: m.V_Tc[t] *( m.Io_Tc[t]/m.a_Tc[t]*pyo.exp( (m.V_Tc[t]+m.I_Tc[t]*m.Rs)/m.a_Tc[t] ) + 1/m.Rsh ) - / ( 1 + m.Rs/m.Rsh + m.Io_Tc[t]*m.Rs/m.a_Tc[t]*pyo.exp( (m.V_Tc[t]+m.I_Tc[t]*m.Rs)/m.a_Tc[t] ) ) - m.I_Tc[t] == 0) -m.f_6 = pyo.Constraint(m.Tc, rule=lambda t: m.Il_Tc[t] - m.Io_Tc[t]*(pyo.exp( (m.V_Tc[t]+m.I_Tc[t]*m.Rs)/m.a_Tc[t] ) - 1) - (m.V_Tc[t] + m.I_Tc[t]*m.Rs)/m.Rsh - m.I_Tc[t]) -m.Pmax_Tc = pyo.Expression(m.Tc, rule=lambda t: m.V_Tc[t] * m.I_Tc[t]) - -def gamma_rule(m, t): - if t >= m.nTc: - return pyo.Constraint.Skip - return (m.Pmax_Tc[t]-m.Pmax_Tc[t-1])*100/(m.Vmp*m.Imp*(m.Tc[t]-m.Tc[t-1])) - -m.gamma_Tc = pyo.Expression(m.Tc, rule=gamma_rule) -m.gamma = pyo.Expression(expr=sum(m.gamma_Tc[t] for t in m.Tc) / m.nTc) - -m.f_7 = pyo.Constraint(expr=m.gamma - m.gPmp == 0) +# Module6ParNonlinear +m.par = pyo.Block() +m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15)) +m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20)) +m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7)) # might need to scale +m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75)) +m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6)) +m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100)) + +m.par.f0 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc == 0) +m.par.f1 = pyo.Constraint(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il == 0) +m.par.f2 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( (m.Vmp + m.Imp*m.par.Rs) / m.par.a ) - 1 ) - (m.Vmp + m.Imp*m.par.Rs)/m.par.Rsh - m.Imp == 0) +m.par.f3 = pyo.Constraint(expr=m.Imp - m.Vmp*( + ( m.par.Io/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + 1/m.par.Rsh ) + /( 1 + m.par.Io*m.par.Rs/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + m.par.Rs/m.par.Rsh ) ) == 0) + +m.par.dT = pyo.Param(initialize=5) + +m.par.aT = pyo.Expression(expr=m.par.a*(m.Tref+m.par.dT)/m.Tref) +m.par.VocT = pyo.Expression(expr=m.bVoc*(1+m.par.Adj/100.0)*m.par.dT + m.Voc) +m.par.Eg = pyo.Expression(expr=(1-0.0002677*m.par.dT)*m.Egref) +m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) +m.par.f4 = pyo.Constraint(expr=m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh == 0) + +# mod6par_gamma_approx +temperatures = np.arange(-10 + 273.15, 50 + 273.15, 3) +def gamma_expr(b, t): + if t == 1: + return pyo.Expression.Skip + return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) + +def gamma_blocks(b, i): + b.Tc = pyo.Param(initialize=temperatures[i - 1]) + b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) + b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) + + # PTnonlinear + b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) + b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) + b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) + b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) + / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) + b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) + b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) + +nTc = len(temperatures) +m.gamma = pyo.Block() +g = m.gamma +g.i = pyo.RangeSet(nTc) +g.pt = pyo.Block(g.i, rule=gamma_blocks) + +g.gamma_Tc = pyo.Expression(g.i, rule=gamma_expr) +g.gamma_avg = pyo.Expression(expr=sum(g.gamma_Tc) / nTc) + +g.f_7 = pyo.Constraint(expr=g.gamma_avg - m.gPmp == 0) # Sanity checks on current +model.sanity = pyo.Block() +s = model.sanity + +s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) +s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) -m.Imp = pyo.Var(domain=pyo.NonNegativeReals) -m.f_8 = pyo.Constraint(expr=m.Il - m.Imp - m.Io * (pyo.exp((m.Vmp + m.Imp * m.Rs) / m.a) - 1.0) - (m.Vmp + m.Imp * m.Rs) / m.Rsh) +s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals) +s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) -m.Ioc = pyo.Var(domain=pyo.NonNegativeReals) -m.f_9 = pyo.Constraint(expr=m.Il - m.Ioc - m.Io * (pyo.exp((m.Voc + m.Ioc * m.Rs) / m.a) - 1.0) - (m.Voc + m.Ioc * m.Rs) / m.Rsh) +# Sanity checks on max_slope +# v_range = np.linspace(0.015 * m.Voc.value, 0.98 * m.Voc.value, 100) -# Sanity checks on slope -v_range = np.linspace(0.015 * m.Voc, 0.98 * m.Voc, 100) +# examine solved modules +model.obj = pyo.Objective(expr=0) +solver = pyo.SolverFactory('ipopt') + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" -with open(filename, 'r') as f: - -if Imp < Isc: - raise ValueError - -Vmp.set_value() -Imp.set_value() -Voc.set_value() -Isc.set_value() -bVoc.set_value() -aIsc.set_value() -gPmp.set_value() -Egref.set_value() -Tref.set_value() +cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + +for n, r in cec_modules_df.iterrows(): + tech = r['Technology'] + + # if r['I_mp_ref'] < r['I_sc_ref']: + # raise ValueError + + m.Vmp.set_value(r['V_mp_ref']) + m.Imp.set_value(r['I_mp_ref']) + m.Voc.set_value(r['V_oc_ref']) + m.Isc.set_value(r['I_sc_ref']) + m.aIsc.set_value(r['alpha_sc']) + m.bVoc.set_value(r['beta_oc']) + m.gPmp.set_value(r['gamma_r']) + m.Tref.set_value(25 + 273.15) + + m.par.a.set_value(r['a_ref']) + m.par.Il.set_value(r['I_L_ref']) + m.par.Io.set_value(r['I_o_ref']) + m.par.Rs.set_value(r['R_s']) + m.par.Rsh.set_value(r['R_sh_ref']) + m.par.Adj.set_value(r['Adjust']) + + log_infeasible_constraints(m.par, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + log_infeasible_bounds(m.par, logger=solve_log, tol=1e-5) + + solver.solve(m.par, tee=True) From 67ba7a99798edac6d820b2b1ab9efe5b7bfe711a Mon Sep 17 00:00:00 2001 From: Darice Date: Wed, 12 Feb 2025 12:54:16 -0700 Subject: [PATCH 03/17] adding stuff --- files/SixParSolve.py | 299 ++++++++++++++++++++++++++++++------------- 1 file changed, 208 insertions(+), 91 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 0dcff7d..ffe8f30 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -6,98 +6,156 @@ solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties") -model = pyo.ConcreteModel() -model.solver = pyo.Block() -m = model.solver - -m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) -m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) -m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) -m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) -m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) -m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - -# Module6ParNonlinear -m.par = pyo.Block() -m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15)) -m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20)) -m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7)) # might need to scale -m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75)) -m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6)) -m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100)) - -m.par.f0 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc == 0) -m.par.f1 = pyo.Constraint(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il == 0) -m.par.f2 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( (m.Vmp + m.Imp*m.par.Rs) / m.par.a ) - 1 ) - (m.Vmp + m.Imp*m.par.Rs)/m.par.Rsh - m.Imp == 0) -m.par.f3 = pyo.Constraint(expr=m.Imp - m.Vmp*( - ( m.par.Io/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + 1/m.par.Rsh ) - /( 1 + m.par.Io*m.par.Rs/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + m.par.Rs/m.par.Rsh ) ) == 0) - -m.par.dT = pyo.Param(initialize=5) - -m.par.aT = pyo.Expression(expr=m.par.a*(m.Tref+m.par.dT)/m.Tref) -m.par.VocT = pyo.Expression(expr=m.bVoc*(1+m.par.Adj/100.0)*m.par.dT + m.Voc) -m.par.Eg = pyo.Expression(expr=(1-0.0002677*m.par.dT)*m.Egref) -m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) -m.par.f4 = pyo.Constraint(expr=m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh == 0) - -# mod6par_gamma_approx -temperatures = np.arange(-10 + 273.15, 50 + 273.15, 3) -def gamma_expr(b, t): - if t == 1: - return pyo.Expression.Skip - return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) - -def gamma_blocks(b, i): - b.Tc = pyo.Param(initialize=temperatures[i - 1]) - b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) - b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) - - # PTnonlinear - b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) - b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) - b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) - b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) - / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) - b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) - b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) - -nTc = len(temperatures) -m.gamma = pyo.Block() -g = m.gamma -g.i = pyo.RangeSet(nTc) -g.pt = pyo.Block(g.i, rule=gamma_blocks) - -g.gamma_Tc = pyo.Expression(g.i, rule=gamma_expr) -g.gamma_avg = pyo.Expression(expr=sum(g.gamma_Tc) / nTc) - -g.f_7 = pyo.Constraint(expr=g.gamma_avg - m.gPmp == 0) - -# Sanity checks on current -model.sanity = pyo.Block() -s = model.sanity - -s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) -s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) - -s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals) -s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) - -# Sanity checks on max_slope -# v_range = np.linspace(0.015 * m.Voc.value, 0.98 * m.Voc.value, 100) - -# examine solved modules -model.obj = pyo.Objective(expr=0) +def plot_iv_curve(model): + curves = [ [ 1000, 25, 'black' ], + [ 800, 25, 'blue' ], + [ 600, 25, 'red' ], + [ 400, 25, 'forest green' ], + [ 200, 25, 'orange' ] ]; + + np = 150 + xx = alloc(np) + yy = alloc(np) + + for( i=0;i<#curves;i++ ) + { + Irr = curves[i][0]; + Tc = curves[i][1]; + cec_model_ivcurve( Irr, Tc, + a, Il, Io, Rs, Rsh, Adj, alpha, I_mp_ref, + vmax, np, + xx, yy ); + + plot(xx,yy,{'series'=sprintf('%g W/m^2', Irr), + size= (Irr==1000 ? 3:1 ), "color"=curves[i][2]}); + + } + + + //Make the plot + plotopt({"title"=sprintf('IV curves at %g ^\\deg C',Tc),"popup"=true,"backcolor"=[255,255,255],"legend"=true, 'legendpos'='bottom'}); + axis('x1', {'label'='Voltage (V)', 'min'=0, 'max'=ceil(vmax*1.01)}); + axis('y1', {'label'='Current (A)', 'min'=0, 'max'= (imax*1.05)}); + +def create_model(): + model = pyo.ConcreteModel() + model.solver = pyo.Block() + m = model.solver + + m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) + m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) + m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) + m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) + m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + + # Module6ParNonlinear + m.par = pyo.Block() + m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15)) + m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20)) + m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7)) # might need to scale + m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75)) + m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6)) + m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100)) + + m.par.f0 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc == 0) + m.par.f1 = pyo.Constraint(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il == 0) + m.par.f2 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( (m.Vmp + m.Imp*m.par.Rs) / m.par.a ) - 1 ) - (m.Vmp + m.Imp*m.par.Rs)/m.par.Rsh - m.Imp == 0) + m.par.f3 = pyo.Constraint(expr=m.Imp - m.Vmp*( + ( m.par.Io/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + 1/m.par.Rsh ) + /( 1 + m.par.Io*m.par.Rs/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + m.par.Rs/m.par.Rsh ) ) == 0) + + m.par.dT = pyo.Param(initialize=5) + + m.par.aT = pyo.Expression(expr=m.par.a*(m.Tref+m.par.dT)/m.Tref) + m.par.VocT = pyo.Expression(expr=m.bVoc*(1+m.par.Adj/100.0)*m.par.dT + m.Voc) + m.par.Eg = pyo.Expression(expr=(1-0.0002677*m.par.dT)*m.Egref) + m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) + m.par.f4 = pyo.Constraint(expr=m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh == 0) + + # mod6par_gamma_approx + temperatures = np.arange(-10 + 273.15, 50 + 273.15, 3) + def gamma_expr(b, t): + if t == 1: + return 0 + return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) + + def gamma_blocks(b, i): + b.Tc = pyo.Param(initialize=temperatures[i - 1]) + b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) + b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) + + # PTnonlinear + b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) + b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) + b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) + b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) + / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) + b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) + b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) + + nTc = len(temperatures) + m.gamma = pyo.Block() + g = m.gamma + g.i = pyo.RangeSet(nTc) + g.pt = pyo.Block(g.i, rule=gamma_blocks) + + g.gamma_Tc = pyo.Expression(g.i, rule=gamma_expr) + g.gamma_avg = pyo.Expression(expr=pyo.summation(g.gamma_Tc) / (nTc - 1)) + + g.f_7 = pyo.Expression(expr=(g.gamma_avg - m.gPmp) ** 2) + + # Sanity checks on current + model.sanity = pyo.Block() + s = model.sanity + + s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) + s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) + + s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals) + s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) + + # Sanity checks on max_slope + # v_range = np.linspace(0.015 * m.Voc.value, 0.98 * m.Voc.value, 100) + model.slope = pyo.Block() + + + # examine solved modules + model.obj = pyo.Objective(rule=g.f_7) + model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + return model solver = pyo.SolverFactory('ipopt') - +solver.options["max_iter"] = 5000 + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) +param_comparison = { + "index": [], + "a_ssc": [], + "Il_ssc": [], + "Io_ssc": [], + "Rs_ssc": [], + "Rsh_ssc": [], + "Adj_ssc": [], + "a_py": [], + "Il_py": [], + "Io_py": [], + "Rs_py": [], + "Rsh_py": [], + "Adj_py": [] +} + +tee = False +IL_SCALING = 1e9 +RSH_SCALING = 1e-3 for n, r in cec_modules_df.iterrows(): + model = create_model() + m = model.solver tech = r['Technology'] # if r['I_mp_ref'] < r['I_sc_ref']: @@ -119,8 +177,67 @@ def gamma_blocks(b, i): m.par.Rsh.set_value(r['R_sh_ref']) m.par.Adj.set_value(r['Adjust']) - log_infeasible_constraints(m.par, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - log_infeasible_bounds(m.par, logger=solve_log, tol=1e-5) - - solver.solve(m.par, tee=True) - + model.scaling_factor[model.solver.par.Io] = IL_SCALING + model.scaling_factor[model.solver.par.Rsh] = RSH_SCALING + scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) + + scaled_model.solver.par.obj = pyo.Objective(rule=0) + + # log_infeasible_constraints(scaled_model.solver.par, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + # log_infeasible_bounds(scaled_model.solver.par, logger=solve_log, tol=1e-5) + + scaled_model.solver.par.del_component('obj') + + solver.solve(scaled_model.solver.par, tee=tee) + + scaled_model.solver.gamma.obj = pyo.Objective(rule=model.solver.gamma.f_7) + + solver.solve(scaled_model.solver.gamma, tee=tee) + # print(pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) + # print(pyo.value(scaled_model.solver.gamma.obj)) + + scaled_model.solver.gamma.del_component('obj') + + solver.solve(scaled_model.sanity, tee=tee) + + try: + res = solver.solve(scaled_model, tee=tee) + except: + pass + + if res.solver.status != 'ok': + # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) + print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) + # print(pyo.value(scaled_model.solver.gamma.obj)) + else: + param_comparison['index'].append(n) + param_comparison['a_ssc'].append(r['a_ref']) + param_comparison['Il_ssc'].append(r['I_L_ref']) + param_comparison['Io_ssc'].append(r['I_o_ref']) + param_comparison['Rs_ssc'].append(r['R_s']) + param_comparison['Rsh_ssc'].append(r['R_sh_ref']) + param_comparison['Adj_ssc'].append(r['Adjust']) + a = pyo.value(scaled_model.solver.par.scaled_a) + Il = pyo.value(scaled_model.solver.par.scaled_Il) + Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(scaled_model.solver.par.scaled_Rs) + Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + + param_comparison['a_py'].append(a) + param_comparison['Il_py'].append(Il) + param_comparison['Io_py'].append(Io) + param_comparison['Rs_py'].append(Rs) + param_comparison['Rsh_py'].append(Rsh) + param_comparison['Adj_py'].append(Adj) + +comparison_df = pd.DataFrame(param_comparison).set_index('index') +comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() +comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() +comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() +comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() +comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() +comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() + +comparison_df.to_csv("param_comparison.csv") \ No newline at end of file From f592ca5611cb08708c3a6969c55592a40610e3ad Mon Sep 17 00:00:00 2001 From: Darice Date: Tue, 18 Mar 2025 13:05:18 -0600 Subject: [PATCH 04/17] add test and plotting --- files/SixParSolve.py | 342 ++++++++++++++++++++++---------------- tests/test_SixParSolve.py | 79 +++++++++ 2 files changed, 280 insertions(+), 141 deletions(-) create mode 100644 tests/test_SixParSolve.py diff --git a/files/SixParSolve.py b/files/SixParSolve.py index ffe8f30..1529eae 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -1,41 +1,99 @@ import pyomo.environ as pyo import pandas as pd import numpy as np +import matplotlib.pyplot as plt import idaes.logger as idaeslog from pyomo.util.infeasible import log_infeasible_constraints, log_infeasible_bounds, log_close_to_bounds solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties") +# 6par_solve.h L 423 +def current_at_voltage_cec(Vmodule, IL_ref, IO_ref, RS, A_ref, RSH_ref, I_mp_ref): + F = 0 + Fprime = 0 + Iold = 0.0 + Inew = I_mp_ref + + it = 0 + maxit = 4000 + while (abs(Inew - Iold) > 1.0e-4 and it < maxit ): + Iold = Inew + + F = IL_ref - Iold - IO_ref * \ + (np.exp((Vmodule + Iold * RS) / A_ref) - 1.0) - \ + (Vmodule + Iold * RS) / RSH_ref + + Fprime = -1.0 - IO_ref * (RS / A_ref) * \ + np.exp((Vmodule + Iold * RS) / A_ref) - \ + (RS / RSH_ref) + + Inew = max(0.0, (Iold - (F / Fprime))) + + return Inew + +def cec_model_params_at_condition(model, Irr, T_cell_K): + Tc_ref = pyo.value(model.Tref) + a = pyo.value(model.par.a) + Io = pyo.value(model.par.Io) + eg0 = pyo.value(model.Egref) + Adj = pyo.value(model.par.Adj) + alpha = pyo.value(model.aIsc) + Il = pyo.value(model.par.Il) + Io = pyo.value(model.par.Io) + Rs = pyo.value(model.par.Rs) + Rsh = pyo.value(model.par.Rsh) + + q = 1.6e-19 + KB = 1.38e-23 + + I_ref = 1000 + muIsc = alpha * (1-Adj/100) + + A_oper = a * T_cell_K / Tc_ref + EG = eg0 * (1-0.0002677*(T_cell_K-Tc_ref)) + # instead of 1/KB, is 11600 in L129 + IO_oper = Io * np.power(T_cell_K/Tc_ref, 3) * np.exp( 11600 *(eg0/Tc_ref - EG/T_cell_K) ) + + Rsh_oper = Rsh*(I_ref/Irr) + IL_oper = Irr/I_ref *( Il + muIsc*(T_cell_K-Tc_ref) ) + if IL_oper < 0.0: + IL_oper = 0.0 + + return IL_oper, IO_oper, Rs, A_oper, Rsh_oper + +def cec_model_ivcurve(model, Irr, T_cell_K, vmax, npts): + I_mp_ref = pyo.value(model.Imp) + IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(model, Irr, T_cell_K) + + y_I = [] + + V = np.linspace(0, vmax, npts) + for i, v in enumerate(V): + I = current_at_voltage_cec( v, IL_oper, IO_oper, Rs, A_oper, Rsh_oper, I_mp_ref ) + y_I.append(I) + return V, y_I + def plot_iv_curve(model): curves = [ [ 1000, 25, 'black' ], - [ 800, 25, 'blue' ], - [ 600, 25, 'red' ], - [ 400, 25, 'forest green' ], - [ 200, 25, 'orange' ] ]; - - np = 150 - xx = alloc(np) - yy = alloc(np) - - for( i=0;i<#curves;i++ ) - { - Irr = curves[i][0]; - Tc = curves[i][1]; - cec_model_ivcurve( Irr, Tc, - a, Il, Io, Rs, Rsh, Adj, alpha, I_mp_ref, - vmax, np, - xx, yy ); - - plot(xx,yy,{'series'=sprintf('%g W/m^2', Irr), - size= (Irr==1000 ? 3:1 ), "color"=curves[i][2]}); - - } - - - //Make the plot - plotopt({"title"=sprintf('IV curves at %g ^\\deg C',Tc),"popup"=true,"backcolor"=[255,255,255],"legend"=true, 'legendpos'='bottom'}); - axis('x1', {'label'='Voltage (V)', 'min'=0, 'max'=ceil(vmax*1.01)}); - axis('y1', {'label'='Current (A)', 'min'=0, 'max'= (imax*1.05)}); + [ 800, 25, 'blue' ], + [ 600, 25, 'red' ], + [ 400, 25, 'green' ], + [ 200, 25, 'orange' ] ] + + alpha = pyo.value(model.aIsc) + beta = pyo.value(model.bVoc) + vmax = pyo.value(model.Voc) + imax = pyo.value(model.Isc) + I_mp_ref = pyo.value(model.Imp) + + npts = 150 + for curve in curves: + Irr = curve[0] + Tc = curve[1] + + x_V, y_I = cec_model_ivcurve(model, Irr, Tc + 273.15, vmax, npts) + plt.plot(x_V, y_I, label=f"{Irr} W/m^2", color=curve[2]) + def create_model(): model = pyo.ConcreteModel() @@ -50,7 +108,7 @@ def create_model(): m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) - m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, units=pyo.units.K) # Module6ParNonlinear m.par = pyo.Block() @@ -128,116 +186,118 @@ def gamma_blocks(b, i): model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) return model -solver = pyo.SolverFactory('ipopt') -solver.options["max_iter"] = 5000 - -filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" -cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) - -param_comparison = { - "index": [], - "a_ssc": [], - "Il_ssc": [], - "Io_ssc": [], - "Rs_ssc": [], - "Rsh_ssc": [], - "Adj_ssc": [], - "a_py": [], - "Il_py": [], - "Io_py": [], - "Rs_py": [], - "Rsh_py": [], - "Adj_py": [] -} - -tee = False -IL_SCALING = 1e9 -RSH_SCALING = 1e-3 -for n, r in cec_modules_df.iterrows(): - model = create_model() - m = model.solver - tech = r['Technology'] - - # if r['I_mp_ref'] < r['I_sc_ref']: - # raise ValueError - - m.Vmp.set_value(r['V_mp_ref']) - m.Imp.set_value(r['I_mp_ref']) - m.Voc.set_value(r['V_oc_ref']) - m.Isc.set_value(r['I_sc_ref']) - m.aIsc.set_value(r['alpha_sc']) - m.bVoc.set_value(r['beta_oc']) - m.gPmp.set_value(r['gamma_r']) - m.Tref.set_value(25 + 273.15) - - m.par.a.set_value(r['a_ref']) - m.par.Il.set_value(r['I_L_ref']) - m.par.Io.set_value(r['I_o_ref']) - m.par.Rs.set_value(r['R_s']) - m.par.Rsh.set_value(r['R_sh_ref']) - m.par.Adj.set_value(r['Adjust']) - - model.scaling_factor[model.solver.par.Io] = IL_SCALING - model.scaling_factor[model.solver.par.Rsh] = RSH_SCALING - scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) - - scaled_model.solver.par.obj = pyo.Objective(rule=0) - - # log_infeasible_constraints(scaled_model.solver.par, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - # log_infeasible_bounds(scaled_model.solver.par, logger=solve_log, tol=1e-5) - - scaled_model.solver.par.del_component('obj') - - solver.solve(scaled_model.solver.par, tee=tee) - - scaled_model.solver.gamma.obj = pyo.Objective(rule=model.solver.gamma.f_7) - - solver.solve(scaled_model.solver.gamma, tee=tee) - # print(pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) - # print(pyo.value(scaled_model.solver.gamma.obj)) - - scaled_model.solver.gamma.del_component('obj') - - solver.solve(scaled_model.sanity, tee=tee) - - try: - res = solver.solve(scaled_model, tee=tee) - except: - pass - - if res.solver.status != 'ok': - # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) - print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) +if __name__ == "__main__": + + solver = pyo.SolverFactory('ipopt') + solver.options["max_iter"] = 5000 + + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" + cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + + param_comparison = { + "index": [], + "a_ssc": [], + "Il_ssc": [], + "Io_ssc": [], + "Rs_ssc": [], + "Rsh_ssc": [], + "Adj_ssc": [], + "a_py": [], + "Il_py": [], + "Io_py": [], + "Rs_py": [], + "Rsh_py": [], + "Adj_py": [] + } + + tee = False + IL_SCALING = 1e9 + RSH_SCALING = 1e-3 + for n, r in cec_modules_df.iterrows(): + model = create_model() + m = model.solver + tech = r['Technology'] + + # if r['I_mp_ref'] < r['I_sc_ref']: + # raise ValueError + + m.Vmp.set_value(r['V_mp_ref']) + m.Imp.set_value(r['I_mp_ref']) + m.Voc.set_value(r['V_oc_ref']) + m.Isc.set_value(r['I_sc_ref']) + m.aIsc.set_value(r['alpha_sc']) + m.bVoc.set_value(r['beta_oc']) + m.gPmp.set_value(r['gamma_r']) + m.Tref.set_value(25 + 273.15) + + m.par.a.set_value(r['a_ref']) + m.par.Il.set_value(r['I_L_ref']) + m.par.Io.set_value(r['I_o_ref']) + m.par.Rs.set_value(r['R_s']) + m.par.Rsh.set_value(r['R_sh_ref']) + m.par.Adj.set_value(r['Adjust']) + + model.scaling_factor[model.solver.par.Io] = IL_SCALING + model.scaling_factor[model.solver.par.Rsh] = RSH_SCALING + scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) + + scaled_model.solver.par.obj = pyo.Objective(rule=0) + + # log_infeasible_constraints(scaled_model.solver.par, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + # log_infeasible_bounds(scaled_model.solver.par, logger=solve_log, tol=1e-5) + + scaled_model.solver.par.del_component('obj') + + solver.solve(scaled_model.solver.par, tee=tee) + + scaled_model.solver.gamma.obj = pyo.Objective(rule=model.solver.gamma.f_7) + + solver.solve(scaled_model.solver.gamma, tee=tee) + # print(pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) # print(pyo.value(scaled_model.solver.gamma.obj)) - else: - param_comparison['index'].append(n) - param_comparison['a_ssc'].append(r['a_ref']) - param_comparison['Il_ssc'].append(r['I_L_ref']) - param_comparison['Io_ssc'].append(r['I_o_ref']) - param_comparison['Rs_ssc'].append(r['R_s']) - param_comparison['Rsh_ssc'].append(r['R_sh_ref']) - param_comparison['Adj_ssc'].append(r['Adjust']) - a = pyo.value(scaled_model.solver.par.scaled_a) - Il = pyo.value(scaled_model.solver.par.scaled_Il) - Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) - Rs = pyo.value(scaled_model.solver.par.scaled_Rs) - Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) - Adj = pyo.value(scaled_model.solver.par.scaled_Adj) - - param_comparison['a_py'].append(a) - param_comparison['Il_py'].append(Il) - param_comparison['Io_py'].append(Io) - param_comparison['Rs_py'].append(Rs) - param_comparison['Rsh_py'].append(Rsh) - param_comparison['Adj_py'].append(Adj) - -comparison_df = pd.DataFrame(param_comparison).set_index('index') -comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() -comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() -comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() -comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() -comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() -comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() - -comparison_df.to_csv("param_comparison.csv") \ No newline at end of file + + scaled_model.solver.gamma.del_component('obj') + + solver.solve(scaled_model.sanity, tee=tee) + + try: + res = solver.solve(scaled_model, tee=tee) + except: + pass + + if res.solver.status != 'ok': + # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) + print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) + # print(pyo.value(scaled_model.solver.gamma.obj)) + else: + param_comparison['index'].append(n) + param_comparison['a_ssc'].append(r['a_ref']) + param_comparison['Il_ssc'].append(r['I_L_ref']) + param_comparison['Io_ssc'].append(r['I_o_ref']) + param_comparison['Rs_ssc'].append(r['R_s']) + param_comparison['Rsh_ssc'].append(r['R_sh_ref']) + param_comparison['Adj_ssc'].append(r['Adjust']) + a = pyo.value(scaled_model.solver.par.scaled_a) + Il = pyo.value(scaled_model.solver.par.scaled_Il) + Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(scaled_model.solver.par.scaled_Rs) + Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + + param_comparison['a_py'].append(a) + param_comparison['Il_py'].append(Il) + param_comparison['Io_py'].append(Io) + param_comparison['Rs_py'].append(Rs) + param_comparison['Rsh_py'].append(Rsh) + param_comparison['Adj_py'].append(Adj) + + comparison_df = pd.DataFrame(param_comparison).set_index('index') + comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() + comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() + comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() + comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() + comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() + comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() + + comparison_df.to_csv("param_comparison.csv") \ No newline at end of file diff --git a/tests/test_SixParSolve.py b/tests/test_SixParSolve.py new file mode 100644 index 0000000..6038977 --- /dev/null +++ b/tests/test_SixParSolve.py @@ -0,0 +1,79 @@ +from pytest import approx +import sys +from pathlib import Path +sys.path.append(str(Path(__file__).parent.parent)) +from files.SixParSolve import * + +def test_default_sam_cec_user(): + # SunPower SPR-E19-310-COM + Vmp = 54.7 + Imp = 5.67 + Voc = 64.4 + Isc = 6.05 + bVoc = -0.175619 + aIsc = 0.00373527 + gamma = -0.386 + # Tref = 46 + 273.15 + Tref = 25 + 273.15 + + a = 2.5776 + Il = 6.054 + Io = 8.360453e-11 + Rs = 3.081202e-01 + Rsh = 500.069 + Adj = 22.909 + + model = create_model() + m = model.solver + m.Vmp.set_value(Vmp) + m.Imp.set_value(Imp) + m.Voc.set_value(Voc) + m.Isc.set_value(Isc) + m.bVoc.set_value(bVoc) + m.aIsc.set_value(aIsc) + m.gPmp.set_value(gamma) + m.Tref.set_value(Tref) + + m.par.a.set_value(a) + m.par.Il.set_value(Il) + m.par.Io.set_value(Io) + m.par.Rs.set_value(Rs) + m.par.Rsh.set_value(Rsh) + m.par.Adj.set_value(Adj) + + IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(m, 1000, 25+275.15) + assert IL_oper == approx(6.059759, rel=1e-5) + assert IO_oper == approx(1.1674203993060455e-10, rel=1e-5) + assert Rs == approx(0.3081202, rel=1e-5) + assert A_oper == approx(2.5948906, rel=1e-5) + assert Rsh_oper == approx(500.069, rel=1e-5) + + plot_iv_curve(m) + + +def test_cec_model_ivcurve_default(): + vmax = 64.4 + muIsc = 0.00287955 + EG = 1.12 + IL_oper = 6.05373 + IO_oper = 8.36045e-11 + Rs = 0.30812 + A_oper = 2.57764 + Rsh_oper = 500.069 + I_mp_ref = 5.67 + + y_I = [] + + V = np.linspace(0, vmax, 150) + for i, v in enumerate(V): + I = current_at_voltage_cec( v, IL_oper, IO_oper, Rs, A_oper, Rsh_oper, I_mp_ref ) + y_I.append(I) + + assert V[1] == approx(0.43221, rel=1e-3) + assert V[-1] == approx(64.4, rel=1e-3) + assert y_I[0] == approx(6.05, rel=1e-3) + assert y_I[1] == approx(6.04913, rel=1e-3) + assert y_I[-1] == approx(5.770225265737295e-06, rel=1e-3) + +test_cec_model_ivcurve_default() +test_default_sam_cec_user() \ No newline at end of file From 2daa4f1ae06ed03d138e622cfc9afa3b39c35795 Mon Sep 17 00:00:00 2001 From: Darice Date: Mon, 21 Apr 2025 09:38:19 -0600 Subject: [PATCH 05/17] solver working for many modules --- files/SixParSolve.py | 496 +++++++++++++++++++++++++++++++++----- tests/test_SixParSolve.py | 61 +++-- 2 files changed, 476 insertions(+), 81 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 1529eae..6bcbe83 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -6,6 +6,9 @@ from pyomo.util.infeasible import log_infeasible_constraints, log_infeasible_bounds, log_close_to_bounds solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties") +IL_SCALING = 1e8 +RSH_SCALING = 1e-3 +param_cols = ["a_py","Il_py","Io_py","Rs_py","Rsh_py","Adj_py"] # 6par_solve.h L 423 def current_at_voltage_cec(Vmodule, IL_ref, IO_ref, RS, A_ref, RSH_ref, I_mp_ref): @@ -32,19 +35,17 @@ def current_at_voltage_cec(Vmodule, IL_ref, IO_ref, RS, A_ref, RSH_ref, I_mp_ref return Inew def cec_model_params_at_condition(model, Irr, T_cell_K): - Tc_ref = pyo.value(model.Tref) - a = pyo.value(model.par.a) - Io = pyo.value(model.par.Io) - eg0 = pyo.value(model.Egref) - Adj = pyo.value(model.par.Adj) - alpha = pyo.value(model.aIsc) - Il = pyo.value(model.par.Il) - Io = pyo.value(model.par.Io) - Rs = pyo.value(model.par.Rs) - Rsh = pyo.value(model.par.Rsh) - - q = 1.6e-19 - KB = 1.38e-23 + m = model.solver + Tc_ref = pyo.value(m.Tref) + a = pyo.value(m.par.a) + Io = pyo.value(m.par.Io) + eg0 = pyo.value(m.Egref) + Adj = pyo.value(m.par.Adj) + alpha = pyo.value(m.aIsc) + Il = pyo.value(m.par.Il) + Io = pyo.value(m.par.Io) + Rs = pyo.value(m.par.Rs) + Rsh = pyo.value(m.par.Rsh) I_ref = 1000 muIsc = alpha * (1-Adj/100) @@ -62,7 +63,7 @@ def cec_model_params_at_condition(model, Irr, T_cell_K): return IL_oper, IO_oper, Rs, A_oper, Rsh_oper def cec_model_ivcurve(model, Irr, T_cell_K, vmax, npts): - I_mp_ref = pyo.value(model.Imp) + I_mp_ref = pyo.value(model.solver.Imp) IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(model, Irr, T_cell_K) y_I = [] @@ -80,11 +81,7 @@ def plot_iv_curve(model): [ 400, 25, 'green' ], [ 200, 25, 'orange' ] ] - alpha = pyo.value(model.aIsc) - beta = pyo.value(model.bVoc) - vmax = pyo.value(model.Voc) - imax = pyo.value(model.Isc) - I_mp_ref = pyo.value(model.Imp) + vmax = pyo.value(model.solver.Voc) npts = 150 for curve in curves: @@ -112,12 +109,12 @@ def create_model(): # Module6ParNonlinear m.par = pyo.Block() - m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15)) - m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20)) - m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7)) # might need to scale - m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75)) - m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6)) - m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100)) + m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15), initialize=7.75) + m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20), initialize=10.25) + m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7), initialize=5e-8) # might need to scale + m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75), initialize=32.5) + m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6), initialize=5e5) + m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100), initialize=1) m.par.f0 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc == 0) m.par.f1 = pyo.Constraint(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il == 0) @@ -182,12 +179,87 @@ def gamma_blocks(b, i): # examine solved modules - model.obj = pyo.Objective(rule=g.f_7) model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) return model -if __name__ == "__main__": +def solve_model(model, solver, tee=False): + model.scaling_factor[model.solver.par.Io] = IL_SCALING + model.scaling_factor[model.solver.par.Rsh] = RSH_SCALING + + scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) + + scaled_model.obj_zero = pyo.Objective(rule=0) + res = None + try: + res = solver.solve(scaled_model, tee=tee) + except Exception as e: + if tee: + log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) + log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) + if res: + return res, scaled_model + else: + return None, scaled_model + if res.solver.status != "ok": + try: + res = solver.solve(scaled_model, tee=tee) + except: + return None, scaled_model + + if res.solver.status != "ok": + scaled_model.solver.gamma.deactivate() + try: + res = solver.solve(scaled_model, tee=tee) + except: + return None, scaled_model + + scaled_model.solver.gamma.activate() + if res.solver.status != "ok": + if tee: + log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) + log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) + if res: + return res, scaled_model + else: + return None, scaled_model + + scaled_model.del_component("obj_zero") + + scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7) + + if pyo.value(scaled_model.solver.gamma.f_7) > abs(pyo.value(scaled_model.solver.gPmp)) * .1: + res = None + try: + res = solver.solve(scaled_model, tee=tee) + except Exception as e: + solver.options["tol"] = 1e-5 + try: + res = solver.solve(scaled_model, tee=tee) + except: + return None, scaled_model + + solver.options["tol"] = 1e-9 + + if tee: + log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) + log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) + + if res.solver.status != "ok": + if tee: + log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) + log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) + if res: + return res, scaled_model + else: + return None, scaled_model + # print("new (gamma - gmp)**2", pyo.value(scaled_model.solver.gamma.f_7), "gPmp", pyo.value(scaled_model.solver.gPmp)) + + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) + return res, scaled_model + + +def run_solved_modules(): solver = pyo.SolverFactory('ipopt') solver.options["max_iter"] = 5000 @@ -211,8 +283,7 @@ def gamma_blocks(b, i): } tee = False - IL_SCALING = 1e9 - RSH_SCALING = 1e-3 + for n, r in cec_modules_df.iterrows(): model = create_model() m = model.solver @@ -221,6 +292,14 @@ def gamma_blocks(b, i): # if r['I_mp_ref'] < r['I_sc_ref']: # raise ValueError + param_comparison['index'].append(n) + param_comparison['a_ssc'].append(r['a_ref']) + param_comparison['Il_ssc'].append(r['I_L_ref']) + param_comparison['Io_ssc'].append(r['I_o_ref']) + param_comparison['Rs_ssc'].append(r['R_s']) + param_comparison['Rsh_ssc'].append(r['R_sh_ref']) + param_comparison['Adj_ssc'].append(r['Adjust']) + m.Vmp.set_value(r['V_mp_ref']) m.Imp.set_value(r['I_mp_ref']) m.Voc.set_value(r['V_oc_ref']) @@ -230,54 +309,264 @@ def gamma_blocks(b, i): m.gPmp.set_value(r['gamma_r']) m.Tref.set_value(25 + 273.15) - m.par.a.set_value(r['a_ref']) - m.par.Il.set_value(r['I_L_ref']) - m.par.Io.set_value(r['I_o_ref']) - m.par.Rs.set_value(r['R_s']) - m.par.Rsh.set_value(r['R_sh_ref']) - m.par.Adj.set_value(r['Adjust']) + # m.par.a.set_value(r['a_ref']) + # m.par.Il.set_value(r['I_L_ref']) + # m.par.Io.set_value(r['I_o_ref']) + # m.par.Rs.set_value(r['R_s']) + # m.par.Rsh.set_value(r['R_sh_ref']) + # m.par.Adj.set_value(r['Adjust']) + + res, scaled_model = solve_model(model, solver, tee) + + if res is None or res.solver.status != 'ok': + # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) + # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) + # print(pyo.value(scaled_model.solver.gamma.obj)) + param_comparison['a_py'].append(None) + param_comparison['Il_py'].append(None) + param_comparison['Io_py'].append(None) + param_comparison['Rs_py'].append(None) + param_comparison['Rsh_py'].append(None) + param_comparison['Adj_py'].append(None) + else: + + a = pyo.value(scaled_model.solver.par.scaled_a) + Il = pyo.value(scaled_model.solver.par.scaled_Il) + Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(scaled_model.solver.par.scaled_Rs) + Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(scaled_model.solver.par.scaled_Adj) - model.scaling_factor[model.solver.par.Io] = IL_SCALING - model.scaling_factor[model.solver.par.Rsh] = RSH_SCALING - scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) + param_comparison['a_py'].append(a) + param_comparison['Il_py'].append(Il) + param_comparison['Io_py'].append(Io) + param_comparison['Rs_py'].append(Rs) + param_comparison['Rsh_py'].append(Rsh) + param_comparison['Adj_py'].append(Adj) - scaled_model.solver.par.obj = pyo.Objective(rule=0) + if n % 10 == 0: + comparison_df = pd.DataFrame(param_comparison).set_index('index') + comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() + comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() + comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() + comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() + comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() + comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() - # log_infeasible_constraints(scaled_model.solver.par, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - # log_infeasible_bounds(scaled_model.solver.par, logger=solve_log, tol=1e-5) + comparison_df.to_csv("param_comparison.csv") - scaled_model.solver.par.del_component('obj') - solver.solve(scaled_model.solver.par, tee=tee) - scaled_model.solver.gamma.obj = pyo.Objective(rule=model.solver.gamma.f_7) +if __name__ == "__main__": - solver.solve(scaled_model.solver.gamma, tee=tee) - # print(pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) - # print(pyo.value(scaled_model.solver.gamma.obj)) + solver = pyo.SolverFactory('ipopt') + solver.options["max_iter"] = 8000 - scaled_model.solver.gamma.del_component('obj') + all_cec_modules_df = pd.read_csv("cec_modules_params.csv") + unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'] == None] - solver.solve(scaled_model.sanity, tee=tee) + df = all_cec_modules_df + model = create_model() + m = model.solver + for i, r in unsolved_df.iterrows(): + diff = ( + (df['V_mp_ref'] - r['V_mp_ref']) ** 2 + + (df['I_mp_ref'] - r['I_mp_ref']) ** 2 + + (df['V_oc_ref'] - r['V_oc_ref']) ** 2 + + (df['I_sc_ref'] - r['I_sc_ref']) ** 2 + + (df['alpha_sc'] - r['alpha_sc']) ** 2 + + (df['beta_oc'] - r['beta_oc']) ** 2 + + (df['gamma_r'] - r['gamma_r']) ** 2) - try: - res = solver.solve(scaled_model, tee=tee) - except: - pass + cec_closest = all_cec_modules_df.iloc[diff.argmin()] + + m.Vmp.set_value(r['V_mp_ref']) + m.Imp.set_value(r['I_mp_ref']) + m.Voc.set_value(r['V_oc_ref']) + m.Isc.set_value(r['I_sc_ref']) + m.aIsc.set_value(r['alpha_sc']) + m.bVoc.set_value(r['beta_oc']) + m.gPmp.set_value(r['gamma_r']) + m.Tref.set_value(25 + 273.15) + + m.par.a.set_value(cec_closest['a_ref']) + m.par.Il.set_value(cec_closest['I_L_ref']) + m.par.Io.set_value(cec_closest['I_o_ref']) + m.par.Rs.set_value(cec_closest['R_s']) + m.par.Rsh.set_value(cec_closest['R_sh_ref']) + m.par.Adj.set_value(cec_closest['Adjust']) + + res, scaled_model = solve_model(model, solver, False) + + if res is None or res.solver.status != 'ok': + continue + else: + a = pyo.value(scaled_model.solver.par.scaled_a) + Il = pyo.value(scaled_model.solver.par.scaled_Il) + Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(scaled_model.solver.par.scaled_Rs) + Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + + all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] + + exit() + + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" + cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/Bad Modules 2024-Nov-15_09-53.csv" + cec_bad_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + cec_bad_modules_df.index = range(len(cec_modules_df), len(cec_modules_df) + len(cec_bad_modules_df)) + + cec_bad_modules_df = cec_bad_modules_df.rename(columns={"Model": "Name", "celltype": "Technology", "vmp": "V_mp_ref", + "imp": "I_mp_ref", "voc": "V_oc_ref", "isc": "I_sc_ref", "alpha_isc": "alpha_sc", + "beta_voc": "beta_oc", "gamma_pmp": "gamma_r", "n_ser": "N_s", "t_ref": "STC"}) + cec_bad_modules_df = cec_bad_modules_df.drop(columns=["Message"]) + + + solved_modules = "/Users/dguittet/Projects/SAM/pysam/param_comparison_1.csv" + solved_modules_df = pd.read_csv(solved_modules, index_col='index').dropna() + cec_modules_solved_df = cec_modules_df[cec_modules_df.index.isin(solved_modules_df.index)] + + solved_bad_modules = "/Users/dguittet/Projects/SAM/pysam/bad_modules_param_0.csv" + solved_bad_modules_df = pd.read_csv(solved_bad_modules, index_col='index') + solved_bad_modules_df.index = range(len(cec_modules_df), len(cec_modules_df) + len(solved_bad_modules_df)) + cec_bad_modules_solved_df = cec_bad_modules_df[cec_bad_modules_df.index.isin(solved_bad_modules_df.index)] + + # read all + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/PV_Module_List_Full_Data_ADA-2024-11-12.xlsx" + all_cec_modules_df = pd.read_excel(filename, skiprows=list(range(0, 16)) + [17]) + all_cec_modules_df = all_cec_modules_df.drop(columns=["Description", 'Safety Certification', + 'Nameplate Pmax', 'Notes', + 'Design Qualification Certification\n(Optional Submission)', + 'Performance Evaluation (Optional Submission)', 'Family', + 'N_p', 'αIpmax', 'βVpmax', 'IPmax, low', 'VPmax, low', 'IPmax, NOCT', + 'VPmax, NOCT', 'Mounting', 'Type', 'Short Side', 'Long Side', + 'Geometric Multiplier', 'P2/Pref', 'CEC Listing Date', 'Last Update']) + all_cec_modules_df = all_cec_modules_df.rename(columns={ + 'Nameplate Isc': "I_sc_ref", 'Nameplate Voc': "V_oc_ref", + 'Nameplate Ipmax': "I_mp_ref", 'Nameplate Vpmax': "V_mp_ref", 'Average NOCT': "T_NOCT", + 'γPmax': "gamma_r", 'αIsc': "alpha_sc", + 'βVoc': "beta_oc", + }) + num_cols = ['A_c', 'N_s', 'I_sc_ref', 'V_oc_ref', 'I_mp_ref', 'V_mp_ref', 'T_NOCT', + 'gamma_r', 'alpha_sc', 'beta_oc'] + all_cec_modules_df[num_cols] = all_cec_modules_df[num_cols].astype(float) + all_cec_modules_df['alpha_sc'] *= all_cec_modules_df['I_sc_ref'] * 1e-2 + all_cec_modules_df['beta_oc'] *= all_cec_modules_df['V_oc_ref'] * 1e-2 + all_cec_modules_df = all_cec_modules_df.drop_duplicates() + + for col in param_cols: + all_cec_modules_df[col] = None + + for i, row in all_cec_modules_df.iterrows(): + df_solved = cec_modules_solved_df[(cec_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) + & (cec_modules_solved_df['Name'].str.contains(row['Model Number'], regex=False)) + & (cec_modules_solved_df['PTC'] == row['PTC'])].drop_duplicates() + if len(df_solved) > 1: + df_solved = df_solved[df_solved['Name'].str.endswith(row['Model Number'])] + if len(df_solved) > 1: + print("Duplicate: ", i, row) + continue + if len(df_solved) == 1: + param_row = solved_modules_df.loc[df_solved.index] + + if len(df_solved) == 0: + df_solved = cec_bad_modules_solved_df[(cec_bad_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) + & (cec_bad_modules_solved_df['Name'].str.contains(row['Model Number'], regex=False)) + ].drop_duplicates() + if len(df_solved) > 1: + df_solved = df_solved[df_solved['Name'].str.endswith(row['Model Number'])] + if len(df_solved) > 1: + print("Duplicate: ", i, row) + continue + if len(df_solved) == 1: + param_row = solved_bad_modules_df.loc[df_solved.index] + + if len(df_solved) == 0: + continue + + all_cec_modules_df.loc[i, param_cols] = param_row[param_cols].values[0] + all_cec_modules_df.to_csv("cec_modules_params.csv") + + exit() + + param_comparison = { + "index": [], + "a_ssc": [], + "Il_ssc": [], + "Io_ssc": [], + "Rs_ssc": [], + "Rsh_ssc": [], + "Adj_ssc": [], + "a_py": [], + "Il_py": [], + "Io_py": [], + "Rs_py": [], + "Rsh_py": [], + "Adj_py": [] + } + + for i, r in cec_modules_unsolved_df.iterrows(): + model = create_model() + m = model.solver + tech = r['Technology'] + + + param_comparison['index'].append(i) + param_comparison['a_ssc'].append(r['a_ref']) + param_comparison['Il_ssc'].append(r['I_L_ref']) + param_comparison['Io_ssc'].append(r['I_o_ref']) + param_comparison['Rs_ssc'].append(r['R_s']) + param_comparison['Rsh_ssc'].append(r['R_sh_ref']) + param_comparison['Adj_ssc'].append(r['Adjust']) + + m.Vmp.set_value(r['V_mp_ref']) + m.Imp.set_value(r['I_mp_ref']) + m.Voc.set_value(r['V_oc_ref']) + m.Isc.set_value(r['I_sc_ref']) + m.aIsc.set_value(r['alpha_sc']) + m.bVoc.set_value(r['beta_oc']) + m.gPmp.set_value(r['gamma_r']) + # m.Tref.set_value(r['STC']) + m.Tref.set_value(25 + 273.15) + + df = cec_modules_solved_df + diff = ( + (df['V_mp_ref'] - r['V_mp_ref']) ** 2 + + (df['I_mp_ref'] - r['I_mp_ref']) ** 2 + + (df['V_oc_ref'] - r['V_oc_ref']) ** 2 + + (df['I_sc_ref'] - r['I_sc_ref']) ** 2 + + (df['alpha_sc'] - r['alpha_sc']) ** 2 + + (df['beta_oc'] - r['beta_oc']) ** 2 + + (df['gamma_r'] - r['gamma_r']) ** 2) + + cec_closest = cec_modules_solved_df.iloc[diff.argmin()] + + m.par.a.set_value(cec_closest['a_ref']) + m.par.Il.set_value(cec_closest['I_L_ref']) + m.par.Io.set_value(cec_closest['I_o_ref']) + m.par.Rs.set_value(cec_closest['R_s']) + m.par.Rsh.set_value(cec_closest['R_sh_ref']) + m.par.Adj.set_value(cec_closest['Adjust']) - if res.solver.status != 'ok': + res, scaled_model = solve_model(model, solver, False) + + if res is None or res.solver.status != 'ok': # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) - print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) + # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) # print(pyo.value(scaled_model.solver.gamma.obj)) + param_comparison['a_py'].append(None) + param_comparison['Il_py'].append(None) + param_comparison['Io_py'].append(None) + param_comparison['Rs_py'].append(None) + param_comparison['Rsh_py'].append(None) + param_comparison['Adj_py'].append(None) else: - param_comparison['index'].append(n) - param_comparison['a_ssc'].append(r['a_ref']) - param_comparison['Il_ssc'].append(r['I_L_ref']) - param_comparison['Io_ssc'].append(r['I_o_ref']) - param_comparison['Rs_ssc'].append(r['R_s']) - param_comparison['Rsh_ssc'].append(r['R_sh_ref']) - param_comparison['Adj_ssc'].append(r['Adjust']) + a = pyo.value(scaled_model.solver.par.scaled_a) Il = pyo.value(scaled_model.solver.par.scaled_Il) Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) @@ -292,6 +581,18 @@ def gamma_blocks(b, i): param_comparison['Rsh_py'].append(Rsh) param_comparison['Adj_py'].append(Adj) + # if i % 100 == 0: + # comparison_df = pd.DataFrame(param_comparison).set_index('index') + # comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() + # comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() + # comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() + # comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() + # comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() + # comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() + + # comparison_df = comparison_df.dropna() + # cec_modules_solved_df = pd.concat([cec_modules_solved_df, comparison_df]) + comparison_df = pd.DataFrame(param_comparison).set_index('index') comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() @@ -300,4 +601,73 @@ def gamma_blocks(b, i): comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() - comparison_df.to_csv("param_comparison.csv") \ No newline at end of file + comparison_df = comparison_df.dropna() + for i, row in comparison_df.iterrows(): + solved_modules_df.loc[i] = row + solved_modules_df.to_csv("param_comparison.csv") + + exit() + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/Bad Modules 2024-Nov-15_09-53.csv" + bad_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + + param_comparison = { + "index": [], + "a_py": [], + "Il_py": [], + "Io_py": [], + "Rs_py": [], + "Rsh_py": [], + "Adj_py": [] + } + + tee = False + + for n, r in bad_modules_df.iterrows(): + model = create_model() + m = model.solver + tech = r['Technology'] + + param_comparison['index'].append(n) + + m.Vmp.set_value(r['vmp']) + m.Imp.set_value(r['imp']) + m.Voc.set_value(r['voc']) + m.Isc.set_value(r['isc']) + m.aIsc.set_value(r['alpha_isc']) + m.bVoc.set_value(r['beta_voc']) + m.gPmp.set_value(r['gamma_pmp']) + m.Tref.set_value(25 + 273.15) + + res, scaled_model = solve_model(model, solver, tee) + + if res is None or res.solver.status != 'ok': + # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) + # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) + # print(pyo.value(scaled_model.solver.gamma.obj)) + param_comparison['a_py'].append(None) + param_comparison['Il_py'].append(None) + param_comparison['Io_py'].append(None) + param_comparison['Rs_py'].append(None) + param_comparison['Rsh_py'].append(None) + param_comparison['Adj_py'].append(None) + else: + + a = pyo.value(scaled_model.solver.par.scaled_a) + Il = pyo.value(scaled_model.solver.par.scaled_Il) + Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(scaled_model.solver.par.scaled_Rs) + Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + + param_comparison['a_py'].append(a) + param_comparison['Il_py'].append(Il) + param_comparison['Io_py'].append(Io) + param_comparison['Rs_py'].append(Rs) + param_comparison['Rsh_py'].append(Rsh) + param_comparison['Adj_py'].append(Adj) + + if n % 10 == 0: + comparison_df = pd.DataFrame(param_comparison).set_index('index') + comparison_df.to_csv("bad_modules_param.csv") + \ No newline at end of file diff --git a/tests/test_SixParSolve.py b/tests/test_SixParSolve.py index 6038977..cbc15cd 100644 --- a/tests/test_SixParSolve.py +++ b/tests/test_SixParSolve.py @@ -4,7 +4,7 @@ sys.path.append(str(Path(__file__).parent.parent)) from files.SixParSolve import * -def test_default_sam_cec_user(): +def define_default_model(set_initial_values=True): # SunPower SPR-E19-310-COM Vmp = 54.7 Imp = 5.67 @@ -16,13 +16,6 @@ def test_default_sam_cec_user(): # Tref = 46 + 273.15 Tref = 25 + 273.15 - a = 2.5776 - Il = 6.054 - Io = 8.360453e-11 - Rs = 3.081202e-01 - Rsh = 500.069 - Adj = 22.909 - model = create_model() m = model.solver m.Vmp.set_value(Vmp) @@ -34,21 +27,32 @@ def test_default_sam_cec_user(): m.gPmp.set_value(gamma) m.Tref.set_value(Tref) - m.par.a.set_value(a) - m.par.Il.set_value(Il) - m.par.Io.set_value(Io) - m.par.Rs.set_value(Rs) - m.par.Rsh.set_value(Rsh) - m.par.Adj.set_value(Adj) + if set_initial_values: + a = 2.5776 + Il = 6.054 + Io = 8.360453e-11 + Rs = 3.081202e-01 + Rsh = 500.069 + Adj = 22.909 + m.par.a.set_value(a) + m.par.Il.set_value(Il) + m.par.Io.set_value(Io) + m.par.Rs.set_value(Rs) + m.par.Rsh.set_value(Rsh) + m.par.Adj.set_value(Adj) + return model - IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(m, 1000, 25+275.15) +def test_default_sam_cec_user(): + model = define_default_model() + + IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(model, 1000, 25+275.15) assert IL_oper == approx(6.059759, rel=1e-5) assert IO_oper == approx(1.1674203993060455e-10, rel=1e-5) assert Rs == approx(0.3081202, rel=1e-5) assert A_oper == approx(2.5948906, rel=1e-5) assert Rsh_oper == approx(500.069, rel=1e-5) - plot_iv_curve(m) + plot_iv_curve(model) def test_cec_model_ivcurve_default(): @@ -75,5 +79,26 @@ def test_cec_model_ivcurve_default(): assert y_I[1] == approx(6.04913, rel=1e-3) assert y_I[-1] == approx(5.770225265737295e-06, rel=1e-3) -test_cec_model_ivcurve_default() -test_default_sam_cec_user() \ No newline at end of file + +def test_cec_model_solve(set_initial_values=True): + model = define_default_model(set_initial_values) + + solver = pyo.SolverFactory('ipopt') + solver.options["max_iter"] = 5000 + solver.options["halt_on_ampl_error"] = 'no' + + res, scaled_model = solve_model(model, solver, tee=True) + a = pyo.value(scaled_model.solver.par.scaled_a) + Il = pyo.value(scaled_model.solver.par.scaled_Il) + Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(scaled_model.solver.par.scaled_Rs) + Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + + plot_iv_curve(model) + +# fig = plt.figure() +# test_default_sam_cec_user() + +# test_cec_model_solve(set_initial_values=False) +# plt.show() \ No newline at end of file From a1366267071c51cd2caf86514ae057c55d0b262f Mon Sep 17 00:00:00 2001 From: Darice Date: Wed, 23 Apr 2025 10:38:58 -0600 Subject: [PATCH 06/17] approx solve is working --- files/SixParSolve.py | 737 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 597 insertions(+), 140 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 6bcbe83..e769b77 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -1,3 +1,5 @@ +from pathlib import Path +import os import pyomo.environ as pyo import pandas as pd import numpy as np @@ -74,22 +76,55 @@ def cec_model_ivcurve(model, Irr, T_cell_K, vmax, npts): y_I.append(I) return V, y_I -def plot_iv_curve(model): +def plot_iv_curve(model, linestyle='solid', label=None, plot_anchors=False): curves = [ [ 1000, 25, 'black' ], [ 800, 25, 'blue' ], [ 600, 25, 'red' ], [ 400, 25, 'green' ], [ 200, 25, 'orange' ] ] + + curves = [ [ 1000, 0, 'black' ], + [ 1000, 25, 'red' ], + [ 1000, 50, 'orange' ], + # [ 800, 25, 'blue' ], + # [ 400, 25, 'green' ], + ] vmax = pyo.value(model.solver.Voc) + alpha_sc = pyo.value(model.solver.aIsc) + I_sc_ref = pyo.value(model.solver.Isc) + beta_oc = pyo.value(model.solver.bVoc) + V_oc_ref = pyo.value(model.solver.Voc) + gamma_r = pyo.value(model.solver.gPmp) + V_mp_ref = pyo.value(model.solver.Vmp) + I_mp_ref = pyo.value(model.solver.Imp) + + # if plot_anchors: + # plt.plot(0, I_sc_ref, marker="o", markersize=15, alpha=0.5) + # plt.plot(V_oc_ref, 0, marker="v", markersize=15, alpha=0.5) + # plt.plot(V_mp_ref, I_mp_ref, marker="*", markersize=15, alpha=0.5) npts = 150 for curve in curves: Irr = curve[0] Tc = curve[1] - x_V, y_I = cec_model_ivcurve(model, Irr, Tc + 273.15, vmax, npts) - plt.plot(x_V, y_I, label=f"{Irr} W/m^2", color=curve[2]) + V_oc = beta_oc * (Tc - 25) + V_oc_ref + x_V, y_I = cec_model_ivcurve(model, Irr, Tc + 273.15, V_oc, npts) + plt.plot(x_V, y_I, label=f"{label} {Irr} W/m^2 {Tc} C", color=curve[2], linestyle=linestyle) + + if plot_anchors and Irr == 1000: + I_sc = alpha_sc * (Tc - 25) + I_sc_ref + plt.plot(0, I_sc, marker="o", markersize=10, alpha=0.5, color=curve[2]) + plt.plot(V_oc, 0, marker="v", markersize=10, alpha=0.5, color=curve[2]) + + P_mp = (gamma_r * (Tc - 25) * 1e-2 + 1) * (V_mp_ref * I_mp_ref) + mp_ind = np.argmax(x_V * y_I) + I_mp = y_I[mp_ind] + V_mp = P_mp / I_mp + # V_mp = x_V[mp_ind] + # I_mp = P_mp / V_mp + plt.plot(V_mp, I_mp, marker="*", markersize=10, alpha=0.5, color=curve[2]) def create_model(): @@ -111,7 +146,7 @@ def create_model(): m.par = pyo.Block() m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15), initialize=7.75) m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20), initialize=10.25) - m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7), initialize=5e-8) # might need to scale + m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1e-13, 1e-7), initialize=5e-8) # might need to scale m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75), initialize=32.5) m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6), initialize=5e5) m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100), initialize=1) @@ -131,11 +166,100 @@ def create_model(): m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) m.par.f4 = pyo.Constraint(expr=m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh == 0) + # mod6par_gamma_approx + temperatures = np.arange(-10 + 273.15, 50 + 273.15, 15) + def gamma_expr(b, t): + return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) + + def gamma_blocks(b, i): + b.Tc = pyo.Param(initialize=temperatures[i - 1]) + b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) + b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) + + # PTnonlinear + b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) + b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) + b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) + b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) + / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) + b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) + b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) + + nTc = len(temperatures) + m.gamma = pyo.Block() + g = m.gamma + g.i = pyo.RangeSet(nTc) + g.d_i = pyo.RangeSet(2, nTc) + g.pt = pyo.Block(g.i, rule=gamma_blocks) + + g.gamma_Tc = pyo.Expression(g.d_i, rule=gamma_expr) + g.gamma_avg = pyo.Expression(expr=pyo.summation(g.gamma_Tc) / len(g.d_i)) + + g.f_7 = pyo.Constraint(expr=(g.gamma_avg - m.gPmp) == 0) + + # Sanity checks on current + model.sanity = pyo.Block() + s = model.sanity + + s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) + s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) + + s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals, initialize=0) + s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) + + # Sanity checks on max_slope + # v_range = np.linspace(0.015 * m.Voc.value, 0.98 * m.Voc.value, 100) + model.slope = pyo.Block() + + + # examine solved modules + model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + return model + + +def create_partly_relaxed_model(): + model = pyo.ConcreteModel() + model.solver = pyo.Block() + m = model.solver + + m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) + m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) + m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) + m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) + m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, units=pyo.units.K) + + # Module6ParNonlinear + m.par = pyo.Block() + m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15), initialize=7.75) + m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20), initialize=10.25) + m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1e-13, 1e-7), initialize=5e-8) # might need to scale + m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75), initialize=32.5) + m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6), initialize=5e5) + m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100), initialize=1) + + m.par.f0 = pyo.Expression(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc) + # m.par.f0 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc == 0) + m.par.f1 = pyo.Constraint(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il == 0) + m.par.f2 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( (m.Vmp + m.Imp*m.par.Rs) / m.par.a ) - 1 ) - (m.Vmp + m.Imp*m.par.Rs)/m.par.Rsh - m.Imp == 0) + m.par.f3 = pyo.Constraint(expr=m.Imp - m.Vmp*( + ( m.par.Io/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + 1/m.par.Rsh ) + /( 1 + m.par.Io*m.par.Rs/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + m.par.Rs/m.par.Rsh ) ) == 0) + + m.par.dT = pyo.Param(initialize=5) + + m.par.aT = pyo.Expression(expr=m.par.a*(m.Tref+m.par.dT)/m.Tref) + m.par.VocT = pyo.Expression(expr=m.bVoc*(1+m.par.Adj/100.0)*m.par.dT + m.Voc) + m.par.Eg = pyo.Expression(expr=(1-0.0002677*m.par.dT)*m.Egref) + m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) + m.par.f4 = pyo.Constraint(expr=m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh == 0) + # mod6par_gamma_approx temperatures = np.arange(-10 + 273.15, 50 + 273.15, 3) def gamma_expr(b, t): - if t == 1: - return 0 return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) def gamma_blocks(b, i): @@ -156,12 +280,15 @@ def gamma_blocks(b, i): m.gamma = pyo.Block() g = m.gamma g.i = pyo.RangeSet(nTc) + g.d_i = pyo.RangeSet(2, nTc) g.pt = pyo.Block(g.i, rule=gamma_blocks) - g.gamma_Tc = pyo.Expression(g.i, rule=gamma_expr) - g.gamma_avg = pyo.Expression(expr=pyo.summation(g.gamma_Tc) / (nTc - 1)) + g.gamma_Tc = pyo.Expression(g.d_i, rule=gamma_expr) + g.gamma_avg = pyo.Expression(expr=pyo.summation(g.gamma_Tc) / len(g.d_i)) g.f_7 = pyo.Expression(expr=(g.gamma_avg - m.gPmp) ** 2) + m.f_0_penalty = pyo.Param(initialize=10, mutable=True) + m.f_0 = pyo.Expression(expr=m.par.f0 ** 2 * m.f_0_penalty) # Sanity checks on current model.sanity = pyo.Block() @@ -170,7 +297,7 @@ def gamma_blocks(b, i): s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) - s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals) + s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals, initialize=0) s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) # Sanity checks on max_slope @@ -182,6 +309,96 @@ def gamma_blocks(b, i): model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) return model +def create_full_relaxed_model(): + model = pyo.ConcreteModel() + model.solver = pyo.Block() + m = model.solver + + m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) + m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) + m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) + m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) + m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) + m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, units=pyo.units.K) + + # Module6ParNonlinear + m.par = pyo.Block() + m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15), initialize=7.75) + m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20), initialize=10.25) + m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7), initialize=5e-8) # might need to scale + m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75), initialize=32.5) + m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6), initialize=5e5) + m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100), initialize=1) + + m.par.f0 = pyo.Expression(expr=(m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc)) + m.par.f1 = pyo.Expression(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il) + m.par.f2 = pyo.Expression(expr=m.par.Il - m.par.Io*( pyo.exp( (m.Vmp + m.Imp*m.par.Rs) / m.par.a ) - 1 ) - (m.Vmp + m.Imp*m.par.Rs)/m.par.Rsh - m.Imp) + m.par.f3 = pyo.Expression(expr=m.Imp - m.Vmp*( + ( m.par.Io/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + 1/m.par.Rsh ) + /( 1 + m.par.Io*m.par.Rs/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + m.par.Rs/m.par.Rsh ) )) + + m.par.dT = pyo.Param(initialize=5) + + m.par.aT = pyo.Expression(expr=m.par.a*(m.Tref+m.par.dT)/m.Tref) + m.par.VocT = pyo.Expression(expr=m.bVoc*(1+m.par.Adj/100.0)*m.par.dT + m.Voc) + m.par.Eg = pyo.Expression(expr=(1-0.0002677*m.par.dT)*m.Egref) + m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) + m.par.f4 = pyo.Expression(expr=(m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh)) + + + # mod6par_gamma_approx + temperatures = np.arange(-10 + 273.15, 50 + 273.15, 3) + def gamma_expr(b, t): + return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) + + def gamma_blocks(b, i): + b.Tc = pyo.Param(initialize=temperatures[i - 1]) + b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) + b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) + + # PTnonlinear + b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) + b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) + b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) + b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) + / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) + b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) + b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) + + nTc = len(temperatures) + m.gamma = pyo.Block() + g = m.gamma + g.i = pyo.RangeSet(nTc) + g.d_i = pyo.RangeSet(2, nTc) + g.pt = pyo.Block(g.i, rule=gamma_blocks) + + g.gamma_Tc = pyo.Expression(g.d_i, rule=gamma_expr) + g.gamma_avg = pyo.Expression(expr=pyo.summation(g.gamma_Tc) / len(g.d_i)) + + g.f_7 = pyo.Expression(expr=(g.gamma_avg - m.gPmp) ** 4) + + m.relaxed_obj = pyo.Objective(expr=(10 * m.par.f0 ** 2 + m.par.f1 ** 2 + m.par.f2 ** 2 + m.par.f3 ** 2 + m.par.f4 ** 2 + m.gamma.f_7) ** 0.5 + # + sum([m.gamma.pt[i].f_5 ** 2 / nTc for i in g.i]) ** 0.5 + # + sum([m.gamma.pt[i].f_6 ** 2 / nTc for i in g.i]) ** 0.5 + ) + + # Sanity checks on current + model.sanity = pyo.Block() + s = model.sanity + + s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) + s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) + + s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals, initialize=0) + s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) + + model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + return model + + def solve_model(model, solver, tee=False): model.scaling_factor[model.solver.par.Io] = IL_SCALING model.scaling_factor[model.solver.par.Rsh] = RSH_SCALING @@ -197,6 +414,7 @@ def solve_model(model, solver, tee=False): log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) if res: + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) return res, scaled_model else: return None, scaled_model @@ -220,13 +438,14 @@ def solve_model(model, solver, tee=False): log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) if res: + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) return res, scaled_model else: return None, scaled_model scaled_model.del_component("obj_zero") - scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7) + scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7 ** 2) if pyo.value(scaled_model.solver.gamma.f_7) > abs(pyo.value(scaled_model.solver.gPmp)) * .1: res = None @@ -250,6 +469,7 @@ def solve_model(model, solver, tee=False): log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) if res: + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) return res, scaled_model else: return None, scaled_model @@ -258,119 +478,159 @@ def solve_model(model, solver, tee=False): pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) return res, scaled_model +def get_iterations(log_file): + with open(log_file, 'r') as f: + for line in f: + if "Number of Iterations....:" in line: + it = line.split(": ")[1] + it_n = int(it) + return it_n + +def get_constraint_infeas(model): + if hasattr(model.solver.par, 'f0'): + vals = (pyo.value(model.solver.par.f0), pyo.value(model.solver.par.f1), pyo.value(model.solver.par.f2), pyo.value(model.solver.par.f3), pyo.value(model.solver.par.f4), pyo.value(model.solver.gamma.f_7)) + else: + if hasattr(model.solver.gamma, 'f_7'): + vals = (pyo.value(model.solver.par.scaled_f0), pyo.value(model.solver.par.scaled_f1), pyo.value(model.solver.par.scaled_f2), pyo.value(model.solver.par.scaled_f3), pyo.value(model.solver.par.scaled_f4), pyo.value(model.solver.gamma.f_7)) + else: + vals = (pyo.value(model.solver.par.scaled_f0), pyo.value(model.solver.par.scaled_f1), pyo.value(model.solver.par.scaled_f2), pyo.value(model.solver.par.scaled_f3), pyo.value(model.solver.par.scaled_f4)) + return [abs(v) for v in vals] -def run_solved_modules(): - solver = pyo.SolverFactory('ipopt') - solver.options["max_iter"] = 5000 +def solve_model_best_solution(model, solver, tee=False): + try: + il_scaling = 10**min(12, -int(np.log10(pyo.value(model.solver.par.Io)))) + rsh_scaling = 10**min(5, -int(np.log10(pyo.value(model.solver.par.Rsh)))) + except: + il_scaling = IL_SCALING + rsh_scaling = RSH_SCALING + model.scaling_factor[model.solver.par.Io] = il_scaling + model.scaling_factor[model.solver.par.Rsh] = rsh_scaling + + scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) - filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" - cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + if hasattr(scaled_model.solver, "f_0"): + scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7 ** 0.5 + scaled_model.solver.f_0 ** 0.5) + elif hasattr(scaled_model.solver.gamma, "f_7"): + scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7 ** 0.5) - param_comparison = { - "index": [], - "a_ssc": [], - "Il_ssc": [], - "Io_ssc": [], - "Rs_ssc": [], - "Rsh_ssc": [], - "Adj_ssc": [], - "a_py": [], - "Il_py": [], - "Io_py": [], - "Rs_py": [], - "Rsh_py": [], - "Adj_py": [] - } + scaled_model.sanity.scaled_f_9.deactivate() + + res = None + try: + res = solver.solve(scaled_model, tee=tee, logfile='ipopt_output.log') + except Exception as e: + pass - tee = False + if res is not None and res.solver.status == "ok": + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) + return res, scaled_model - for n, r in cec_modules_df.iterrows(): - model = create_model() - m = model.solver - tech = r['Technology'] + it_n = get_iterations('ipopt_output.log') - # if r['I_mp_ref'] < r['I_sc_ref']: - # raise ValueError + solver.options['max_iter'] = it_n - 1 - param_comparison['index'].append(n) - param_comparison['a_ssc'].append(r['a_ref']) - param_comparison['Il_ssc'].append(r['I_L_ref']) - param_comparison['Io_ssc'].append(r['I_o_ref']) - param_comparison['Rs_ssc'].append(r['R_s']) - param_comparison['Rsh_ssc'].append(r['R_sh_ref']) - param_comparison['Adj_ssc'].append(r['Adjust']) + try: + res = solver.solve(scaled_model, tee=tee, logfile='ipopt_output.log') + except: + pass - m.Vmp.set_value(r['V_mp_ref']) - m.Imp.set_value(r['I_mp_ref']) - m.Voc.set_value(r['V_oc_ref']) - m.Isc.set_value(r['I_sc_ref']) - m.aIsc.set_value(r['alpha_sc']) - m.bVoc.set_value(r['beta_oc']) - m.gPmp.set_value(r['gamma_r']) - m.Tref.set_value(25 + 273.15) + while 'infeasible' in res.solver.message: + it_n = get_iterations('ipopt_output.log') + solver.options['max_iter'] = it_n - 1 + res = solver.solve(scaled_model, tee=tee, logfile='ipopt_output.log') - # m.par.a.set_value(r['a_ref']) - # m.par.Il.set_value(r['I_L_ref']) - # m.par.Io.set_value(r['I_o_ref']) - # m.par.Rs.set_value(r['R_s']) - # m.par.Rsh.set_value(r['R_sh_ref']) - # m.par.Adj.set_value(r['Adjust']) + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - res, scaled_model = solve_model(model, solver, tee) + # get somewhat stable params + infeas_old = np.inf - if res is None or res.solver.status != 'ok': - # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) - # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) - # print(pyo.value(scaled_model.solver.gamma.obj)) - param_comparison['a_py'].append(None) - param_comparison['Il_py'].append(None) - param_comparison['Io_py'].append(None) - param_comparison['Rs_py'].append(None) - param_comparison['Rsh_py'].append(None) - param_comparison['Adj_py'].append(None) - else: + infeas = sum(get_constraint_infeas(model)) + # while infeas_old >= infeas: + attempts = 0 + while infeas > 10 and attempts < 10: + res = solver.solve(scaled_model, tee=tee, logfile='ipopt_output.log') + infeas = sum(get_constraint_infeas(scaled_model)) + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - a = pyo.value(scaled_model.solver.par.scaled_a) - Il = pyo.value(scaled_model.solver.par.scaled_Il) - Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) - Rs = pyo.value(scaled_model.solver.par.scaled_Rs) - Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) - Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + return res, scaled_model - param_comparison['a_py'].append(a) - param_comparison['Il_py'].append(Il) - param_comparison['Io_py'].append(Io) - param_comparison['Rs_py'].append(Rs) - param_comparison['Rsh_py'].append(Rsh) - param_comparison['Adj_py'].append(Adj) - if n % 10 == 0: - comparison_df = pd.DataFrame(param_comparison).set_index('index') - comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() - comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() - comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() - comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() - comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() - comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() +def solve_relaxed_model(model, solver, tee=False): + model.scaling_factor[model.solver.par.Io] = 10**min(12, -int(np.log10(pyo.value(model.solver.par.Io)))) + model.scaling_factor[model.solver.par.Rsh] = 10**min(5, -int(np.log10(pyo.value(model.solver.par.Rsh)))) + + scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) - comparison_df.to_csv("param_comparison.csv") + res = None + solver.options['max_iter'] = 3000 + try: + res = solver.solve(scaled_model, tee=tee) + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) + + except Exception as e: + if tee: + log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) + log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) + if res: + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) + return res, scaled_model + else: + return None, scaled_model + if res.solver.status != "ok": + try: + res = solver.solve(scaled_model, tee=tee) + except: + return None, scaled_model + if res.solver.status != "ok": + scaled_model.solver.gamma.deactivate() + try: + res = solver.solve(scaled_model, tee=tee) + except: + return None, scaled_model -if __name__ == "__main__": + scaled_model.solver.gamma.activate() + if res.solver.status != "ok": + if tee: + log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) + log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) + if res: + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) + return res, scaled_model + else: + return None, scaled_model - solver = pyo.SolverFactory('ipopt') - solver.options["max_iter"] = 8000 + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) + return res, scaled_model - all_cec_modules_df = pd.read_csv("cec_modules_params.csv") - unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'] == None] +def set_parameters(m, r): + try: + m.Vmp.set_value(r['V_mp_ref']) + m.Imp.set_value(r['I_mp_ref']) + m.Voc.set_value(r['V_oc_ref']) + m.Isc.set_value(r['I_sc_ref']) + m.aIsc.set_value(r['alpha_sc']) + m.bVoc.set_value(r['beta_oc']) + m.gPmp.set_value(r['gamma_r']) + m.Tref.set_value(25 + 273.15) + return True + except: + return False - df = all_cec_modules_df - model = create_model() +def set_initial_guess(model, a, Il, Io, Rs, Rsh, Adj): m = model.solver - for i, r in unsolved_df.iterrows(): - diff = ( + m.par.a.set_value(a) + m.par.Il.set_value(Il) + m.par.Io.set_value(Io) + m.par.Rs.set_value(Rs) + m.par.Rsh.set_value(Rsh) + m.par.Adj.set_value(Adj) + model.sanity.Imp_calc.set_value(pyo.value(m.Imp)) + # model.sanity.Ioc_calc.set_value(pyo.value(m.par.c)) + +def find_closest(df, r): + diff = ( (df['V_mp_ref'] - r['V_mp_ref']) ** 2 + (df['I_mp_ref'] - r['I_mp_ref']) ** 2 + (df['V_oc_ref'] - r['V_oc_ref']) ** 2 + @@ -379,38 +639,138 @@ def run_solved_modules(): (df['beta_oc'] - r['beta_oc']) ** 2 + (df['gamma_r'] - r['gamma_r']) ** 2) - cec_closest = all_cec_modules_df.iloc[diff.argmin()] + cec_closest = df.iloc[diff.argmin()] + return cec_closest + +def get_params_from_model(model): + if hasattr(model.solver.par, 'scaled_a'): + a = pyo.value(model.solver.par.scaled_a) + Il = pyo.value(model.solver.par.scaled_Il) + Io = pyo.value(model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(model.solver.par.scaled_Rs) + Rsh = pyo.value(model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(model.solver.par.scaled_Adj) + else: + a = pyo.value(model.solver.par.a) + Il = pyo.value(model.solver.par.Il) + Io = pyo.value(model.solver.par.Io) + Rs = pyo.value(model.solver.par.Rs) + Rsh = pyo.value(model.solver.par.Rsh) + Adj = pyo.value(model.solver.par.Adj) + return a, Il, Io, Rs, Rsh, Adj + + +def copy_over_vars(model, scaled_model, model_rel): + pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - m.Vmp.set_value(r['V_mp_ref']) - m.Imp.set_value(r['I_mp_ref']) - m.Voc.set_value(r['V_oc_ref']) - m.Isc.set_value(r['I_sc_ref']) - m.aIsc.set_value(r['alpha_sc']) - m.bVoc.set_value(r['beta_oc']) - m.gPmp.set_value(r['gamma_r']) - m.Tref.set_value(25 + 273.15) + set_initial_guess(model_rel, *get_params_from_model(scaled_model)) + for i, gamma_block in model.solver.gamma.pt.items(): + model_rel.solver.gamma.pt[i].V_Tc.set_value(pyo.value(gamma_block.V_Tc)) + model_rel.solver.gamma.pt[i].I_Tc.set_value(pyo.value(gamma_block.I_Tc)) - m.par.a.set_value(cec_closest['a_ref']) - m.par.Il.set_value(cec_closest['I_L_ref']) - m.par.Io.set_value(cec_closest['I_o_ref']) - m.par.Rs.set_value(cec_closest['R_s']) - m.par.Rsh.set_value(cec_closest['R_sh_ref']) - m.par.Adj.set_value(cec_closest['Adjust']) +def get_curve_diffs(r, model): + x_V, y_I = cec_model_ivcurve(model, 1000, 25 + 273.15, r['V_oc_ref'], 150) + p = x_V * y_I + mp_ind = np.argmax(p) + d_I_sc = (y_I[0] - r['I_sc_ref']) / r['I_sc_ref'] + d_I_mp = (y_I[mp_ind] - r['I_mp_ref']) / r['I_mp_ref'] + d_V_mp = (x_V[mp_ind] - r['V_mp_ref']) / r['V_mp_ref'] + d_P_mp = (p[mp_ind] - r['V_mp_ref'] * r['I_mp_ref']) / (r['V_mp_ref'] * r['I_mp_ref']) + return d_I_sc, d_I_mp, d_V_mp, d_P_mp - res, scaled_model = solve_model(model, solver, False) - if res is None or res.solver.status != 'ok': +if __name__ == "__main__": + + output_path = Path(__file__).parent / "6parsolve_output" + if not (output_path).exists(): + os.mkdir(output_path) + + solver = pyo.SolverFactory('ipopt') + + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" + sam_solved_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + + all_cec_modules_df = pd.read_csv("cec_modules_params_approx.csv") + unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] + + solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] + + diff_cols = ['d_Isc', 'd_Imp', 'd_Vmp', 'd_Pmp'] + for i, r in unsolved_df.iterrows(): + if r['V_mp_ref'] < 0: + all_cec_modules_df['Error'] = "Vmp < 0" + continue + if r['V_oc_ref'] < r['V_mp_ref']: + all_cec_modules_df['Error'] = "Voc < Vmp" continue + + model = create_model() + if not set_parameters(model.solver, r): + all_cec_modules_df['Error'] = "Parameter missing or out of bounds" + continue + + plt.figure() + solver.options["max_iter"] = 3000 + + # plot SAM one + r_sam = sam_solved_modules_df[(sam_solved_modules_df['Manufacturer'] == r['Manufacturer'].replace(",", "")) + & (sam_solved_modules_df['Name'].str.contains(r['Model Number'], regex=False)) + & (sam_solved_modules_df['PTC'] == r['PTC'])].drop_duplicates() + if len(r_sam): + r_sam = r_sam.to_dict('records')[0] + set_initial_guess(model, r_sam['a_ref'], r_sam['I_L_ref'], r_sam['I_o_ref'], r_sam['R_s'], r_sam['R_sh_ref'], r_sam['Adjust']) + plot_iv_curve(model, linestyle='dashed', label="SAM") + if True: + cec_closest = find_closest(solved_df, r) + set_initial_guess(model, *cec_closest[param_cols]) + plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") + + res, scaled_model = solve_model_best_solution(model, solver, True) + + print(get_constraint_infeas(model)) + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) else: - a = pyo.value(scaled_model.solver.par.scaled_a) - Il = pyo.value(scaled_model.solver.par.scaled_Il) - Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) - Rs = pyo.value(scaled_model.solver.par.scaled_Rs) - Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) - Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) + + # if res is None or res.solver.status != 'ok': + # solver.options['max_iter'] = 1938 + # res = solver.solve(scaled_model, tee=True) + + # model_rel = create_full_relaxed_model() + # set_parameters(model_rel.solver, r) + + # copy_over_vars(model, scaled_model, model_rel) + + # res_rel, scaled_model_rel = solve_relaxed_model(model_rel, solver, True) + + # print(get_constraint_infeas(model_rel)) + # plot_iv_curve(model_rel, linestyle='-', label="Closest", plot_anchors=True) + + # if res_rel is None or res_rel.solver.status != 'ok': + # continue + + d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + + infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() + infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + + if infeas_max > 0.5 or infeas_sum > 0.5: + continue + + a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) + all_cec_modules_df.loc[i, diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] + all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] - all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(output_path / f"IV_curve_{'nosam' if not len(r_sam) else '_'}{i}.png") + plt.close() + all_cec_modules_df.to_csv("cec_modules_params_approx.csv") exit() filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" @@ -462,30 +822,30 @@ def run_solved_modules(): all_cec_modules_df[col] = None for i, row in all_cec_modules_df.iterrows(): - df_solved = cec_modules_solved_df[(cec_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) + r_sam = cec_modules_solved_df[(cec_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) & (cec_modules_solved_df['Name'].str.contains(row['Model Number'], regex=False)) & (cec_modules_solved_df['PTC'] == row['PTC'])].drop_duplicates() - if len(df_solved) > 1: - df_solved = df_solved[df_solved['Name'].str.endswith(row['Model Number'])] - if len(df_solved) > 1: + if len(r_sam) > 1: + r_sam = r_sam[r_sam['Name'].str.endswith(row['Model Number'])] + if len(r_sam) > 1: print("Duplicate: ", i, row) continue - if len(df_solved) == 1: - param_row = solved_modules_df.loc[df_solved.index] + if len(r_sam) == 1: + param_row = solved_modules_df.loc[r_sam.index] - if len(df_solved) == 0: - df_solved = cec_bad_modules_solved_df[(cec_bad_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) + if len(r_sam) == 0: + r_sam = cec_bad_modules_solved_df[(cec_bad_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) & (cec_bad_modules_solved_df['Name'].str.contains(row['Model Number'], regex=False)) ].drop_duplicates() - if len(df_solved) > 1: - df_solved = df_solved[df_solved['Name'].str.endswith(row['Model Number'])] - if len(df_solved) > 1: + if len(r_sam) > 1: + r_sam = r_sam[r_sam['Name'].str.endswith(row['Model Number'])] + if len(r_sam) > 1: print("Duplicate: ", i, row) continue - if len(df_solved) == 1: - param_row = solved_bad_modules_df.loc[df_solved.index] + if len(r_sam) == 1: + param_row = solved_bad_modules_df.loc[r_sam.index] - if len(df_solved) == 0: + if len(r_sam) == 0: continue all_cec_modules_df.loc[i, param_cols] = param_row[param_cols].values[0] @@ -670,4 +1030,101 @@ def run_solved_modules(): if n % 10 == 0: comparison_df = pd.DataFrame(param_comparison).set_index('index') comparison_df.to_csv("bad_modules_param.csv") - \ No newline at end of file + + +def run_solved_modules(): + solver = pyo.SolverFactory('ipopt') + solver.options["max_iter"] = 5000 + + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" + cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) + + param_comparison = { + "index": [], + "a_ssc": [], + "Il_ssc": [], + "Io_ssc": [], + "Rs_ssc": [], + "Rsh_ssc": [], + "Adj_ssc": [], + "a_py": [], + "Il_py": [], + "Io_py": [], + "Rs_py": [], + "Rsh_py": [], + "Adj_py": [] + } + + tee = False + + for n, r in cec_modules_df.iterrows(): + model = create_model() + m = model.solver + tech = r['Technology'] + + # if r['I_mp_ref'] < r['I_sc_ref']: + # raise ValueError + + param_comparison['index'].append(n) + param_comparison['a_ssc'].append(r['a_ref']) + param_comparison['Il_ssc'].append(r['I_L_ref']) + param_comparison['Io_ssc'].append(r['I_o_ref']) + param_comparison['Rs_ssc'].append(r['R_s']) + param_comparison['Rsh_ssc'].append(r['R_sh_ref']) + param_comparison['Adj_ssc'].append(r['Adjust']) + + m.Vmp.set_value(r['V_mp_ref']) + m.Imp.set_value(r['I_mp_ref']) + m.Voc.set_value(r['V_oc_ref']) + m.Isc.set_value(r['I_sc_ref']) + m.aIsc.set_value(r['alpha_sc']) + m.bVoc.set_value(r['beta_oc']) + m.gPmp.set_value(r['gamma_r']) + m.Tref.set_value(25 + 273.15) + + # m.par.a.set_value(r['a_ref']) + # m.par.Il.set_value(r['I_L_ref']) + # m.par.Io.set_value(r['I_o_ref']) + # m.par.Rs.set_value(r['R_s']) + # m.par.Rsh.set_value(r['R_sh_ref']) + # m.par.Adj.set_value(r['Adjust']) + + res, scaled_model = solve_model(model, solver, tee) + + if res is None or res.solver.status != 'ok': + # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) + # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) + # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) + # print(pyo.value(scaled_model.solver.gamma.obj)) + param_comparison['a_py'].append(None) + param_comparison['Il_py'].append(None) + param_comparison['Io_py'].append(None) + param_comparison['Rs_py'].append(None) + param_comparison['Rsh_py'].append(None) + param_comparison['Adj_py'].append(None) + else: + + a = pyo.value(scaled_model.solver.par.scaled_a) + Il = pyo.value(scaled_model.solver.par.scaled_Il) + Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) + Rs = pyo.value(scaled_model.solver.par.scaled_Rs) + Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) + Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + + param_comparison['a_py'].append(a) + param_comparison['Il_py'].append(Il) + param_comparison['Io_py'].append(Io) + param_comparison['Rs_py'].append(Rs) + param_comparison['Rsh_py'].append(Rsh) + param_comparison['Adj_py'].append(Adj) + + if n % 10 == 0: + comparison_df = pd.DataFrame(param_comparison).set_index('index') + comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() + comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() + comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() + comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() + comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() + comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() + + comparison_df.to_csv("param_comparison.csv") From 9116b61b02fd27709ff3c91b122382cdd602e1a9 Mon Sep 17 00:00:00 2001 From: Darice Date: Wed, 30 Apr 2025 13:52:30 -0600 Subject: [PATCH 07/17] cec solver working for most modules --- files/SixParSolve.py | 208 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 175 insertions(+), 33 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index e769b77..022c3b1 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -83,11 +83,11 @@ def plot_iv_curve(model, linestyle='solid', label=None, plot_anchors=False): [ 400, 25, 'green' ], [ 200, 25, 'orange' ] ] - curves = [ [ 1000, 0, 'black' ], + curves = [ [ 1000, 2, 'black' ], [ 1000, 25, 'red' ], - [ 1000, 50, 'orange' ], - # [ 800, 25, 'blue' ], - # [ 400, 25, 'green' ], + [ 1000, 47, 'orange' ], + [ 800, 0, 'blue' ], + [ 400, 0, 'green' ], ] vmax = pyo.value(model.solver.Voc) @@ -127,7 +127,7 @@ def plot_iv_curve(model, linestyle='solid', label=None, plot_anchors=False): plt.plot(V_mp, I_mp, marker="*", markersize=10, alpha=0.5, color=curve[2]) -def create_model(): +def create_model(gamma_curve_dt=3): model = pyo.ConcreteModel() model.solver = pyo.Block() m = model.solver @@ -167,7 +167,7 @@ def create_model(): m.par.f4 = pyo.Constraint(expr=m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh == 0) # mod6par_gamma_approx - temperatures = np.arange(-10 + 273.15, 50 + 273.15, 15) + temperatures = np.arange(-10 + 273.15, 50 + 273.15, gamma_curve_dt) def gamma_expr(b, t): return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) @@ -443,11 +443,11 @@ def solve_model(model, solver, tee=False): else: return None, scaled_model - scaled_model.del_component("obj_zero") + # scaled_model.del_component("obj_zero") - scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7 ** 2) + # scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7 ** 2) - if pyo.value(scaled_model.solver.gamma.f_7) > abs(pyo.value(scaled_model.solver.gPmp)) * .1: + if res is None or res.solver.status != 'ok': res = None try: res = solver.solve(scaled_model, tee=tee) @@ -679,6 +679,127 @@ def get_curve_diffs(r, model): return d_I_sc, d_I_mp, d_V_mp, d_P_mp +def format_cec_datasheet(filename): + all_cec_modules_df = pd.read_excel(filename, skiprows=list(range(0, 16)) + [17]) + all_cec_modules_df = all_cec_modules_df.drop(columns=["Description", 'Safety Certification', + 'Nameplate Pmax', 'Notes', + 'Design Qualification Certification\n(Optional Submission)', + 'Performance Evaluation (Optional Submission)', 'Family', + 'N_p', 'αIpmax', 'βVpmax', 'IPmax, low', 'VPmax, low', 'IPmax, NOCT', + 'VPmax, NOCT', 'Mounting', 'Type', 'Short Side', 'Long Side', + 'Geometric Multiplier', 'P2/Pref', 'CEC Listing Date', 'Last Update']) + all_cec_modules_df = all_cec_modules_df.rename(columns={ + 'Nameplate Isc': "I_sc_ref", 'Nameplate Voc': "V_oc_ref", + 'Nameplate Ipmax': "I_mp_ref", 'Nameplate Vpmax': "V_mp_ref", 'Average NOCT': "T_NOCT", + 'γPmax': "gamma_r", 'αIsc': "alpha_sc", + 'βVoc': "beta_oc", + }) + num_cols = ['A_c', 'N_s', 'I_sc_ref', 'V_oc_ref', 'I_mp_ref', 'V_mp_ref', 'T_NOCT', + 'gamma_r', 'alpha_sc', 'beta_oc'] + all_cec_modules_df[num_cols] = all_cec_modules_df[num_cols].astype(float) + all_cec_modules_df['alpha_sc'] *= all_cec_modules_df['I_sc_ref'] * 1e-2 + all_cec_modules_df['beta_oc'] *= all_cec_modules_df['V_oc_ref'] * 1e-2 + all_cec_modules_df = all_cec_modules_df.drop_duplicates() + return all_cec_modules_df + + +def solve_all_without_guess(all_cec_modules_df, plotting): + for i, r in all_cec_modules_df.iterrows(): + if r['V_mp_ref'] < 0: + all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" + continue + if r['V_oc_ref'] < r['V_mp_ref']: + all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" + continue + + model = create_model() + if not set_parameters(model.solver, r): + all_cec_modules_df.loc[i, 'Error'] = "Parameter missing or out of bounds" + continue + + if plotting: + plt.figure() + solver.options["max_iter"] = 3000 + + res, scaled_model = solve_model(model, solver, tee=False) + + print(get_constraint_infeas(model)) + if plotting: + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) + else: + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) + + if res is None or res.solver.status != 'ok': + continue + + a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) + all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] + all_cec_modules_df.to_csv("cec_modules_params.csv", index=False) + + +def solve_all_with_guess(all_cec_modules_df, plotting): + diff_cols = ['d_Isc', 'd_Imp', 'd_Vmp', 'd_Pmp'] + + unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] + solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] + for i, r in unsolved_df.iterrows(): + if r['V_mp_ref'] < 0: + all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" + continue + if r['V_oc_ref'] < r['V_mp_ref']: + all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" + continue + + model = create_model(gamma_curve_dt=15) + if not set_parameters(model.solver, r): + all_cec_modules_df.loc[i, 'Error'] = "Parameter missing or out of bounds" + continue + + if plotting: + plt.figure() + solver.options["max_iter"] = 3000 + + cec_closest = find_closest(solved_df, r) + set_initial_guess(model, *cec_closest[param_cols]) + if plotting: + plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") + + res, scaled_model = solve_model_best_solution(model, solver, tee=False) + + print(get_constraint_infeas(model)) + if plotting: + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) + else: + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) + + d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + + infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() + infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + + if infeas_max > 0.5 or infeas_sum > 0.5: + all_cec_modules_df.loc[i, 'Error'] = "Infeasibility > 0.5" + if plotting: + plt.close() + continue + + a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) + all_cec_modules_df.loc[i, diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] + all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] + + if plotting: + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(output_path / f"IV_curve_{'nosam' if not len(r_sam) else '_'}{i}.png") + plt.close() + + all_cec_modules_df.to_csv("cec_modules_params_approx.csv", index=False) + if __name__ == "__main__": output_path = Path(__file__).parent / "6parsolve_output" @@ -691,25 +812,29 @@ def get_curve_diffs(r, model): sam_solved_modules_df = pd.read_csv(filename, skiprows=[1, 2]) all_cec_modules_df = pd.read_csv("cec_modules_params_approx.csv") + unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] diff_cols = ['d_Isc', 'd_Imp', 'd_Vmp', 'd_Pmp'] - for i, r in unsolved_df.iterrows(): + + plotting = True + for i, r in all_cec_modules_df.iterrows(): if r['V_mp_ref'] < 0: - all_cec_modules_df['Error'] = "Vmp < 0" + all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" continue if r['V_oc_ref'] < r['V_mp_ref']: - all_cec_modules_df['Error'] = "Voc < Vmp" + all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" continue - model = create_model() + model = create_model(gamma_curve_dt=15) if not set_parameters(model.solver, r): - all_cec_modules_df['Error'] = "Parameter missing or out of bounds" + all_cec_modules_df.loc[i, 'Error'] = "Parameter missing or out of bounds" continue - plt.figure() + if plotting: + plt.figure() solver.options["max_iter"] = 3000 # plot SAM one @@ -719,19 +844,23 @@ def get_curve_diffs(r, model): if len(r_sam): r_sam = r_sam.to_dict('records')[0] set_initial_guess(model, r_sam['a_ref'], r_sam['I_L_ref'], r_sam['I_o_ref'], r_sam['R_s'], r_sam['R_sh_ref'], r_sam['Adjust']) - plot_iv_curve(model, linestyle='dashed', label="SAM") + if plotting: + plot_iv_curve(model, linestyle='dashed', label="SAM") if True: cec_closest = find_closest(solved_df, r) set_initial_guess(model, *cec_closest[param_cols]) - plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") + if plotting: + plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") - res, scaled_model = solve_model_best_solution(model, solver, True) + # res, scaled_model = solve_model_best_solution(model, solver, tee=False) + res, scaled_model = solve_model(model, solver, tee=False) print(get_constraint_infeas(model)) - if res and res.solver.status == 'ok': - plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) - else: - plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) + if plotting: + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) + else: + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) # if res is None or res.solver.status != 'ok': # solver.options['max_iter'] = 1938 @@ -756,23 +885,32 @@ def get_curve_diffs(r, model): infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() if infeas_max > 0.5 or infeas_sum > 0.5: + all_cec_modules_df.loc[i, 'Error'] = "Infeasibility > 0.5" + if plotting: + plt.close() continue a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) all_cec_modules_df.loc[i, diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] - plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - plt.xlabel("Voltage") - plt.ylabel("Current") - plt.title(f"{r['Manufacturer']} {r['Model Number']}") - plt.tight_layout() - plt.savefig(output_path / f"IV_curve_{'nosam' if not len(r_sam) else '_'}{i}.png") - plt.close() + if plotting: + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(output_path / f"IV_curve_{'nosam' if not len(r_sam) else '_'}{i}.png") + plt.close() - all_cec_modules_df.to_csv("cec_modules_params_approx.csv") + all_cec_modules_df.to_csv("cec_modules_params_approx.csv", index=False) exit() + + + + + filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) @@ -853,6 +991,10 @@ def get_curve_diffs(r, model): exit() + + + + param_comparison = { "index": [], "a_ssc": [], @@ -912,7 +1054,7 @@ def get_curve_diffs(r, model): m.par.Rsh.set_value(cec_closest['R_sh_ref']) m.par.Adj.set_value(cec_closest['Adjust']) - res, scaled_model = solve_model(model, solver, False) + res, scaled_model = solve_model(model, solver, tee=False) if res is None or res.solver.status != 'ok': # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) @@ -998,7 +1140,7 @@ def get_curve_diffs(r, model): m.gPmp.set_value(r['gamma_pmp']) m.Tref.set_value(25 + 273.15) - res, scaled_model = solve_model(model, solver, tee) + res, scaled_model = solve_model(model, solver, tee=tee) if res is None or res.solver.status != 'ok': # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) @@ -1089,7 +1231,7 @@ def run_solved_modules(): # m.par.Rsh.set_value(r['R_sh_ref']) # m.par.Adj.set_value(r['Adjust']) - res, scaled_model = solve_model(model, solver, tee) + res, scaled_model = solve_model(model, solver, tee=tee) if res is None or res.solver.status != 'ok': # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) From 0d7ec83515d91fe6e9be143b7d16e0d10faabe68 Mon Sep 17 00:00:00 2001 From: Darice Date: Tue, 27 May 2025 11:51:18 -0600 Subject: [PATCH 08/17] finalize SixParSolve.py as run script --- files/SixParSolve.py | 1008 +++++++++++++----------------------------- 1 file changed, 305 insertions(+), 703 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 022c3b1..3af4d26 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -1,19 +1,30 @@ from pathlib import Path import os +import sys import pyomo.environ as pyo import pandas as pd import numpy as np import matplotlib.pyplot as plt import idaes.logger as idaeslog from pyomo.util.infeasible import log_infeasible_constraints, log_infeasible_bounds, log_close_to_bounds +import logging +from datetime import datetime +logging.getLogger('pyomo.core').setLevel(logging.ERROR) solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties") IL_SCALING = 1e8 RSH_SCALING = 1e-3 -param_cols = ["a_py","Il_py","Io_py","Rs_py","Rsh_py","Adj_py"] +test_data_cols = ['A_c', 'N_s', 'I_sc_ref', 'V_oc_ref', 'I_mp_ref', 'V_mp_ref', 'T_NOCT', + 'gamma_r', 'alpha_sc', 'beta_oc'] +model_param_cols = ["a_py", "Il_py", "Io_py", "Rs_py", "Rsh_py", "Adj_py"] +model_param_cols_ssc = ["a_ssc", "Il_ssc", "Io_ssc","Rs_ssc","Rsh_ssc","Adj_ssc"] +iv_diff_cols = ['d_Isc', 'd_Imp', 'd_Vmp', 'd_Pmp'] # 6par_solve.h L 423 def current_at_voltage_cec(Vmodule, IL_ref, IO_ref, RS, A_ref, RSH_ref, I_mp_ref): + """ + Solve for the current at the voltage Vmodule, given the single diode parameters which could have been modified for non-STC condition + """ F = 0 Fprime = 0 Iold = 0.0 @@ -37,6 +48,9 @@ def current_at_voltage_cec(Vmodule, IL_ref, IO_ref, RS, A_ref, RSH_ref, I_mp_ref return Inew def cec_model_params_at_condition(model, Irr, T_cell_K): + """ + Get the single diode parameters at a given condition + """ m = model.solver Tc_ref = pyo.value(m.Tref) a = pyo.value(m.par.a) @@ -64,7 +78,11 @@ def cec_model_params_at_condition(model, Irr, T_cell_K): return IL_oper, IO_oper, Rs, A_oper, Rsh_oper + def cec_model_ivcurve(model, Irr, T_cell_K, vmax, npts): + """ + Calculate the IV curve with voltage on x-axis and current on y-axis at the given irradiance and temperature + """ I_mp_ref = pyo.value(model.solver.Imp) IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(model, Irr, T_cell_K) @@ -77,6 +95,11 @@ def cec_model_ivcurve(model, Irr, T_cell_K, vmax, npts): return V, y_I def plot_iv_curve(model, linestyle='solid', label=None, plot_anchors=False): + """ + Plot IV curves for a range of conditions using a model from `create_model` with all the test and model parameters populated + + Model must have populated data from `test_data_cols` and `model_param_cols` + """ curves = [ [ 1000, 25, 'black' ], [ 800, 25, 'blue' ], [ 600, 25, 'red' ], @@ -99,11 +122,6 @@ def plot_iv_curve(model, linestyle='solid', label=None, plot_anchors=False): V_mp_ref = pyo.value(model.solver.Vmp) I_mp_ref = pyo.value(model.solver.Imp) - # if plot_anchors: - # plt.plot(0, I_sc_ref, marker="o", markersize=15, alpha=0.5) - # plt.plot(V_oc_ref, 0, marker="v", markersize=15, alpha=0.5) - # plt.plot(V_mp_ref, I_mp_ref, marker="*", markersize=15, alpha=0.5) - npts = 150 for curve in curves: Irr = curve[0] @@ -120,14 +138,21 @@ def plot_iv_curve(model, linestyle='solid', label=None, plot_anchors=False): P_mp = (gamma_r * (Tc - 25) * 1e-2 + 1) * (V_mp_ref * I_mp_ref) mp_ind = np.argmax(x_V * y_I) - I_mp = y_I[mp_ind] - V_mp = P_mp / I_mp - # V_mp = x_V[mp_ind] - # I_mp = P_mp / V_mp + if Tc == 25: + I_mp = I_mp_ref + V_mp = V_mp_ref + else: + I_mp = y_I[mp_ind] + V_mp = P_mp / I_mp plt.plot(V_mp, I_mp, marker="*", markersize=10, alpha=0.5, color=curve[2]) def create_model(gamma_curve_dt=3): + """ + Create a pyomo model for a single diode model with a set of non-STC temperatures. + The set of non-STC temperatures are used to fit the model to the gamma Pmp test parameter + This set of temperatures is sampled using the `gamma_curve_dt` parameter, which is the interval between each temperature sample + """ model = pyo.ConcreteModel() model.solver = pyo.Block() m = model.solver @@ -136,16 +161,16 @@ def create_model(gamma_curve_dt=3): m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) - m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) - m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) + m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True, units=pyo.units.A/pyo.units.K) + m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True, units=pyo.units.V/pyo.units.K) + m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True, units=pyo.units.percent/pyo.units.K) m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, units=pyo.units.K) # Module6ParNonlinear m.par = pyo.Block() m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15), initialize=7.75) - m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20), initialize=10.25) + m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.01, 20), initialize=10.25) m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1e-13, 1e-7), initialize=5e-8) # might need to scale m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75), initialize=32.5) m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6), initialize=5e5) @@ -169,21 +194,22 @@ def create_model(gamma_curve_dt=3): # mod6par_gamma_approx temperatures = np.arange(-10 + 273.15, 50 + 273.15, gamma_curve_dt) def gamma_expr(b, t): - return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) + return (b.pt[t].Pmp_Tc-b.pt[t-1].Pmp_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) def gamma_blocks(b, i): b.Tc = pyo.Param(initialize=temperatures[i - 1]) - b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) - b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) + b.Vmp_Tc = pyo.Var(domain=pyo.NonNegativeReals) + b.Imp_Tc = pyo.Var(domain=pyo.NonNegativeReals) # PTnonlinear b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) - b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) + b.Eg_Tc = pyo.Expression(expr=m.Egref * (1-0.0002677*(b.Tc-m.Tref))) + b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - b.Eg_Tc/b.Tc))) b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) - b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) - / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) - b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) - b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) + b.f_5 = pyo.Constraint(expr=b.Vmp_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.Vmp_Tc+b.Imp_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) + / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.Vmp_Tc+b.Imp_Tc*m.par.Rs)/b.a_Tc ) ) - b.Imp_Tc == 0) + b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.Vmp_Tc+b.Imp_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.Vmp_Tc + b.Imp_Tc*m.par.Rs)/m.par.Rsh - b.Imp_Tc == 0) + b.Pmp_Tc = pyo.Expression(expr=b.Vmp_Tc * b.Imp_Tc) nTc = len(temperatures) m.gamma = pyo.Block() @@ -211,195 +237,19 @@ def gamma_blocks(b, i): # v_range = np.linspace(0.015 * m.Voc.value, 0.98 * m.Voc.value, 100) model.slope = pyo.Block() - - # examine solved modules - model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) - return model - - -def create_partly_relaxed_model(): - model = pyo.ConcreteModel() - model.solver = pyo.Block() - m = model.solver - - m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) - m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) - m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) - m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) - m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, units=pyo.units.K) - - # Module6ParNonlinear - m.par = pyo.Block() - m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15), initialize=7.75) - m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20), initialize=10.25) - m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1e-13, 1e-7), initialize=5e-8) # might need to scale - m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75), initialize=32.5) - m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6), initialize=5e5) - m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100), initialize=1) - - m.par.f0 = pyo.Expression(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc) - # m.par.f0 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc == 0) - m.par.f1 = pyo.Constraint(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il == 0) - m.par.f2 = pyo.Constraint(expr=m.par.Il - m.par.Io*( pyo.exp( (m.Vmp + m.Imp*m.par.Rs) / m.par.a ) - 1 ) - (m.Vmp + m.Imp*m.par.Rs)/m.par.Rsh - m.Imp == 0) - m.par.f3 = pyo.Constraint(expr=m.Imp - m.Vmp*( - ( m.par.Io/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + 1/m.par.Rsh ) - /( 1 + m.par.Io*m.par.Rs/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + m.par.Rs/m.par.Rsh ) ) == 0) - - m.par.dT = pyo.Param(initialize=5) - - m.par.aT = pyo.Expression(expr=m.par.a*(m.Tref+m.par.dT)/m.Tref) - m.par.VocT = pyo.Expression(expr=m.bVoc*(1+m.par.Adj/100.0)*m.par.dT + m.Voc) - m.par.Eg = pyo.Expression(expr=(1-0.0002677*m.par.dT)*m.Egref) - m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) - m.par.f4 = pyo.Constraint(expr=m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh == 0) - - # mod6par_gamma_approx - temperatures = np.arange(-10 + 273.15, 50 + 273.15, 3) - def gamma_expr(b, t): - return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) - - def gamma_blocks(b, i): - b.Tc = pyo.Param(initialize=temperatures[i - 1]) - b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) - b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) - - # PTnonlinear - b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) - b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) - b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) - b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) - / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) - b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) - b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) - - nTc = len(temperatures) - m.gamma = pyo.Block() - g = m.gamma - g.i = pyo.RangeSet(nTc) - g.d_i = pyo.RangeSet(2, nTc) - g.pt = pyo.Block(g.i, rule=gamma_blocks) - - g.gamma_Tc = pyo.Expression(g.d_i, rule=gamma_expr) - g.gamma_avg = pyo.Expression(expr=pyo.summation(g.gamma_Tc) / len(g.d_i)) - - g.f_7 = pyo.Expression(expr=(g.gamma_avg - m.gPmp) ** 2) - m.f_0_penalty = pyo.Param(initialize=10, mutable=True) - m.f_0 = pyo.Expression(expr=m.par.f0 ** 2 * m.f_0_penalty) - - # Sanity checks on current - model.sanity = pyo.Block() - s = model.sanity - - s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) - s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) - - s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals, initialize=0) - s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) - - # Sanity checks on max_slope - # v_range = np.linspace(0.015 * m.Voc.value, 0.98 * m.Voc.value, 100) - model.slope = pyo.Block() - - # examine solved modules model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) return model -def create_full_relaxed_model(): - model = pyo.ConcreteModel() - model.solver = pyo.Block() - m = model.solver - - m.Vmp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.Imp = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.Voc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.Isc = pyo.Param(domain=pyo.NonNegativeReals, mutable=True) - m.aIsc = pyo.Param(domain=pyo.Reals, mutable=True) - m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True) - m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True) - m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) - m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, units=pyo.units.K) - - # Module6ParNonlinear - m.par = pyo.Block() - m.par.a = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.05, 15), initialize=7.75) - m.par.Il = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.5, 20), initialize=10.25) - m.par.Io = pyo.Var(domain=pyo.NonNegativeReals, bounds=(None, 1e-7), initialize=5e-8) # might need to scale - m.par.Rs = pyo.Var(domain=pyo.NonNegativeReals, bounds=(0.001, 75), initialize=32.5) - m.par.Rsh = pyo.Var(domain=pyo.NonNegativeReals, bounds=(1, 1e6), initialize=5e5) - m.par.Adj = pyo.Var(domain=pyo.Reals, bounds=(-100, 100), initialize=1) - - m.par.f0 = pyo.Expression(expr=(m.par.Il - m.par.Io*( pyo.exp( m.Isc*m.par.Rs / m.par.a ) - 1 ) - m.Isc*m.par.Rs/m.par.Rsh - m.Isc)) - m.par.f1 = pyo.Expression(expr=m.par.Io*( pyo.exp( m.Voc/m.par.a ) - 1 ) + m.Voc/m.par.Rsh -m.par.Il) - m.par.f2 = pyo.Expression(expr=m.par.Il - m.par.Io*( pyo.exp( (m.Vmp + m.Imp*m.par.Rs) / m.par.a ) - 1 ) - (m.Vmp + m.Imp*m.par.Rs)/m.par.Rsh - m.Imp) - m.par.f3 = pyo.Expression(expr=m.Imp - m.Vmp*( - ( m.par.Io/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + 1/m.par.Rsh ) - /( 1 + m.par.Io*m.par.Rs/m.par.a*pyo.exp( (m.Vmp + m.Imp*m.par.Rs)/m.par.a ) + m.par.Rs/m.par.Rsh ) )) - - m.par.dT = pyo.Param(initialize=5) - - m.par.aT = pyo.Expression(expr=m.par.a*(m.Tref+m.par.dT)/m.Tref) - m.par.VocT = pyo.Expression(expr=m.bVoc*(1+m.par.Adj/100.0)*m.par.dT + m.Voc) - m.par.Eg = pyo.Expression(expr=(1-0.0002677*m.par.dT)*m.Egref) - m.par.IoT = pyo.Expression(expr=m.par.Io*( (m.Tref+m.par.dT)/m.Tref )**3 *pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/(m.Tref+m.par.dT)))) - m.par.f4 = pyo.Expression(expr=(m.par.Il+m.aIsc*(1-m.par.Adj/100)*m.par.dT - m.par.IoT*(pyo.exp( m.par.VocT/m.par.aT ) - 1 ) - m.par.VocT/m.par.Rsh)) - - - # mod6par_gamma_approx - temperatures = np.arange(-10 + 273.15, 50 + 273.15, 3) - def gamma_expr(b, t): - return (b.pt[t].Pmax_Tc-b.pt[t-1].Pmax_Tc)*100/(m.Vmp*m.Imp*(b.pt[t].Tc-b.pt[t-1].Tc)) - - def gamma_blocks(b, i): - b.Tc = pyo.Param(initialize=temperatures[i - 1]) - b.V_Tc = pyo.Var(domain=pyo.NonNegativeReals) - b.I_Tc = pyo.Var(domain=pyo.NonNegativeReals) - - # PTnonlinear - b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) - b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - m.par.Eg/b.Tc))) - b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) - b.f_5 = pyo.Constraint(expr=b.V_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) - / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) ) - b.I_Tc == 0) - b.f_6 = pyo.Constraint(expr=b.Il_Tc - b.Io_Tc*(pyo.exp( (b.V_Tc+b.I_Tc*m.par.Rs)/b.a_Tc ) - 1) - (b.V_Tc + b.I_Tc*m.par.Rs)/m.par.Rsh - b.I_Tc == 0) - b.Pmax_Tc = pyo.Expression(expr=b.V_Tc * b.I_Tc) - - nTc = len(temperatures) - m.gamma = pyo.Block() - g = m.gamma - g.i = pyo.RangeSet(nTc) - g.d_i = pyo.RangeSet(2, nTc) - g.pt = pyo.Block(g.i, rule=gamma_blocks) - - g.gamma_Tc = pyo.Expression(g.d_i, rule=gamma_expr) - g.gamma_avg = pyo.Expression(expr=pyo.summation(g.gamma_Tc) / len(g.d_i)) - - g.f_7 = pyo.Expression(expr=(g.gamma_avg - m.gPmp) ** 4) - - m.relaxed_obj = pyo.Objective(expr=(10 * m.par.f0 ** 2 + m.par.f1 ** 2 + m.par.f2 ** 2 + m.par.f3 ** 2 + m.par.f4 ** 2 + m.gamma.f_7) ** 0.5 - # + sum([m.gamma.pt[i].f_5 ** 2 / nTc for i in g.i]) ** 0.5 - # + sum([m.gamma.pt[i].f_6 ** 2 / nTc for i in g.i]) ** 0.5 - ) - - # Sanity checks on current - model.sanity = pyo.Block() - s = model.sanity - - s.Imp_calc = pyo.Var(domain=pyo.NonNegativeReals) - s.f_8 = pyo.Constraint(expr=m.par.Il - s.Imp_calc - m.par.Io * (pyo.exp((m.Vmp + s.Imp_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Vmp + s.Imp_calc * m.par.Rs) / m.par.Rsh == 0) - - s.Ioc_calc = pyo.Var(domain=pyo.NonNegativeReals, initialize=0) - s.f_9 = pyo.Constraint(expr=m.par.Il - s.Ioc_calc - m.par.Io * (pyo.exp((m.Voc + s.Ioc_calc * m.par.Rs) / m.par.a) - 1.0) - (m.Voc + s.Ioc_calc * m.par.Rs) / m.par.Rsh == 0) - - model.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) - return model - def solve_model(model, solver, tee=False): + """ + Solve the model with scaling factors, multiple tries and separating steps + + Solution may not be optimal! Solutions may have slight infeasibility. + Caller needs to check whether it is above an acceptable threshold using the log functions by setting tee=True, + or after the function, using `get_constraint_infeas` or `get_curve_diffs`. + """ model.scaling_factor[model.solver.par.Io] = IL_SCALING model.scaling_factor[model.solver.par.Rsh] = RSH_SCALING @@ -443,10 +293,6 @@ def solve_model(model, solver, tee=False): else: return None, scaled_model - # scaled_model.del_component("obj_zero") - - # scaled_model.obj_gamma = pyo.Objective(rule=scaled_model.solver.gamma.f_7 ** 2) - if res is None or res.solver.status != 'ok': res = None try: @@ -473,12 +319,15 @@ def solve_model(model, solver, tee=False): return res, scaled_model else: return None, scaled_model - # print("new (gamma - gmp)**2", pyo.value(scaled_model.solver.gamma.f_7), "gPmp", pyo.value(scaled_model.solver.gPmp)) pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) return res, scaled_model + def get_iterations(log_file): + """ + Get the number of IPOPT iterations from the log file + """ with open(log_file, 'r') as f: for line in f: if "Number of Iterations....:" in line: @@ -486,7 +335,11 @@ def get_iterations(log_file): it_n = int(it) return it_n + def get_constraint_infeas(model): + """ + Get the magnitude of infeasibility for each constraint in the pyomo model. Model can be the original or the scaled version + """ if hasattr(model.solver.par, 'f0'): vals = (pyo.value(model.solver.par.f0), pyo.value(model.solver.par.f1), pyo.value(model.solver.par.f2), pyo.value(model.solver.par.f3), pyo.value(model.solver.par.f4), pyo.value(model.solver.gamma.f_7)) else: @@ -496,7 +349,15 @@ def get_constraint_infeas(model): vals = (pyo.value(model.solver.par.scaled_f0), pyo.value(model.solver.par.scaled_f1), pyo.value(model.solver.par.scaled_f2), pyo.value(model.solver.par.scaled_f3), pyo.value(model.solver.par.scaled_f4)) return [abs(v) for v in vals] + def solve_model_best_solution(model, solver, tee=False): + """ + Solve the model and return the best solution regardless of regardless of whether IPOPT has converged or exited gracefully + + Solution may not be optimal! Solutions may have a lot of infeasibility. + Caller needs to check whether it is above an acceptable threshold using the log functions by setting tee=True, + or after the function, using `get_constraint_infeas` or `get_curve_diffs`. + """ try: il_scaling = 10**min(12, -int(np.log10(pyo.value(model.solver.par.Io)))) rsh_scaling = 10**min(5, -int(np.log10(pyo.value(model.solver.par.Rsh)))) @@ -542,69 +403,21 @@ def solve_model_best_solution(model, solver, tee=False): pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) # get somewhat stable params - infeas_old = np.inf - infeas = sum(get_constraint_infeas(model)) - # while infeas_old >= infeas: attempts = 0 while infeas > 10 and attempts < 10: res = solver.solve(scaled_model, tee=tee, logfile='ipopt_output.log') infeas = sum(get_constraint_infeas(scaled_model)) pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) + attempts += 1 return res, scaled_model -def solve_relaxed_model(model, solver, tee=False): - model.scaling_factor[model.solver.par.Io] = 10**min(12, -int(np.log10(pyo.value(model.solver.par.Io)))) - model.scaling_factor[model.solver.par.Rsh] = 10**min(5, -int(np.log10(pyo.value(model.solver.par.Rsh)))) - - scaled_model = pyo.TransformationFactory('core.scale_model').create_using(model) - - res = None - solver.options['max_iter'] = 3000 - try: - res = solver.solve(scaled_model, tee=tee) - pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - - except Exception as e: - if tee: - log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) - log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) - if res: - pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - return res, scaled_model - else: - return None, scaled_model - - if res.solver.status != "ok": - try: - res = solver.solve(scaled_model, tee=tee) - except: - return None, scaled_model - - if res.solver.status != "ok": - scaled_model.solver.gamma.deactivate() - try: - res = solver.solve(scaled_model, tee=tee) - except: - return None, scaled_model - - scaled_model.solver.gamma.activate() - if res.solver.status != "ok": - if tee: - log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-7) - log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-7) - if res: - pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - return res, scaled_model - else: - return None, scaled_model - - pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - return res, scaled_model - -def set_parameters(m, r): +def set_parameters(m, r: pd.Series): + """ + Set the test data parameters for the model.solver block from a pandas Series + """ try: m.Vmp.set_value(r['V_mp_ref']) m.Imp.set_value(r['I_mp_ref']) @@ -619,6 +432,9 @@ def set_parameters(m, r): return False def set_initial_guess(model, a, Il, Io, Rs, Rsh, Adj): + """ + Set the model parameters from a guess + """ m = model.solver m.par.a.set_value(a) m.par.Il.set_value(Il) @@ -627,22 +443,57 @@ def set_initial_guess(model, a, Il, Io, Rs, Rsh, Adj): m.par.Rsh.set_value(Rsh) m.par.Adj.set_value(Adj) model.sanity.Imp_calc.set_value(pyo.value(m.Imp)) - # model.sanity.Ioc_calc.set_value(pyo.value(m.par.c)) -def find_closest(df, r): +# intercept, coefficients +Voc_Vmp_to_a = [0.22328699, 0.03106774, 0.00738466] +Isc_Imp_to_Il = [0.0153132, 0.96058469, 0.04133798] +Isc_Imp_Voc_Vmp_to_Rs = [0.28969838, 0.61785177, -0.6759943, 0.02241771, 0.0218622] + +def set_empirical_initial_guess(model): + """ + Set the model parameters from an empirical equation based on previously-solved modules + """ + m = model.solver + + Vmp = pyo.value(m.Vmp) + Imp = pyo.value(m.Imp) + Voc = pyo.value(m.Voc) + Isc = pyo.value(m.Isc) + aIsc = pyo.value(m.aIsc) + bVoc = pyo.value(m.bVoc) + gPmp = pyo.value(m.gPmp) + + a = Voc_Vmp_to_a[0] + Voc_Vmp_to_a[1] * Voc + Voc_Vmp_to_a[2] * Vmp + Il = Isc_Imp_to_Il[0] + Isc_Imp_to_Il[1] * Isc + Isc_Imp_to_Il[2] * Imp + Rs = Isc_Imp_Voc_Vmp_to_Rs[0] + Isc_Imp_Voc_Vmp_to_Rs[1] * Isc + Isc_Imp_Voc_Vmp_to_Rs[2] * Imp \ + + Isc_Imp_Voc_Vmp_to_Rs[3] * Voc + Isc_Imp_Voc_Vmp_to_Rs[4] * Vmp + + m.par.a.set_value(a) + m.par.Il.set_value(Il) + m.par.Rs.set_value(Rs) + + +def find_closest(df_solved: pd.DataFrame, r: pd.Series): + """ + Given a dataframe of solved modules, find the module that is closest to the one in row r and return the parameters + """ diff = ( - (df['V_mp_ref'] - r['V_mp_ref']) ** 2 + - (df['I_mp_ref'] - r['I_mp_ref']) ** 2 + - (df['V_oc_ref'] - r['V_oc_ref']) ** 2 + - (df['I_sc_ref'] - r['I_sc_ref']) ** 2 + - (df['alpha_sc'] - r['alpha_sc']) ** 2 + - (df['beta_oc'] - r['beta_oc']) ** 2 + - (df['gamma_r'] - r['gamma_r']) ** 2) + (df_solved['V_mp_ref'] - r['V_mp_ref']) ** 2 + + (df_solved['I_mp_ref'] - r['I_mp_ref']) ** 2 + + (df_solved['V_oc_ref'] - r['V_oc_ref']) ** 2 + + (df_solved['I_sc_ref'] - r['I_sc_ref']) ** 2 + + (df_solved['alpha_sc'] - r['alpha_sc']) ** 2 + + (df_solved['beta_oc'] - r['beta_oc']) ** 2 + + (df_solved['gamma_r'] - r['gamma_r']) ** 2) - cec_closest = df.iloc[diff.argmin()] - return cec_closest + params_closest = df_solved.iloc[diff.argmin()] + return params_closest + def get_params_from_model(model): + """ + Extract the model parameters from the pyomo model + """ if hasattr(model.solver.par, 'scaled_a'): a = pyo.value(model.solver.par.scaled_a) Il = pyo.value(model.solver.par.scaled_Il) @@ -660,15 +511,10 @@ def get_params_from_model(model): return a, Il, Io, Rs, Rsh, Adj -def copy_over_vars(model, scaled_model, model_rel): - pyo.TransformationFactory('core.scale_model').propagate_solution(scaled_model, model) - - set_initial_guess(model_rel, *get_params_from_model(scaled_model)) - for i, gamma_block in model.solver.gamma.pt.items(): - model_rel.solver.gamma.pt[i].V_Tc.set_value(pyo.value(gamma_block.V_Tc)) - model_rel.solver.gamma.pt[i].I_Tc.set_value(pyo.value(gamma_block.I_Tc)) - def get_curve_diffs(r, model): + """ + Calculate the normalized difference in I_sc, I_mp, V_mp and P_mp at 1000 irradiance and 25 C + """ x_V, y_I = cec_model_ivcurve(model, 1000, 25 + 273.15, r['V_oc_ref'], 150) p = x_V * y_I mp_ind = np.argmax(p) @@ -679,8 +525,11 @@ def get_curve_diffs(r, model): return d_I_sc, d_I_mp, d_V_mp, d_P_mp -def format_cec_datasheet(filename): - all_cec_modules_df = pd.read_excel(filename, skiprows=list(range(0, 16)) + [17]) +def read_prepare_file(xlsx_file): + """ + Read the CEC Module Excel Spreadsheet and prepare the file + """ + all_cec_modules_df = pd.read_excel(xlsx_file, skiprows=list(range(0, 16)) + [17]) all_cec_modules_df = all_cec_modules_df.drop(columns=["Description", 'Safety Certification', 'Nameplate Pmax', 'Notes', 'Design Qualification Certification\n(Optional Submission)', @@ -694,16 +543,22 @@ def format_cec_datasheet(filename): 'γPmax': "gamma_r", 'αIsc': "alpha_sc", 'βVoc': "beta_oc", }) - num_cols = ['A_c', 'N_s', 'I_sc_ref', 'V_oc_ref', 'I_mp_ref', 'V_mp_ref', 'T_NOCT', - 'gamma_r', 'alpha_sc', 'beta_oc'] - all_cec_modules_df[num_cols] = all_cec_modules_df[num_cols].astype(float) + all_cec_modules_df[test_data_cols] = all_cec_modules_df[test_data_cols].astype(float) all_cec_modules_df['alpha_sc'] *= all_cec_modules_df['I_sc_ref'] * 1e-2 all_cec_modules_df['beta_oc'] *= all_cec_modules_df['V_oc_ref'] * 1e-2 all_cec_modules_df = all_cec_modules_df.drop_duplicates() + + for col in model_param_cols: + all_cec_modules_df[col] = None return all_cec_modules_df -def solve_all_without_guess(all_cec_modules_df, plotting): +def run_solve_first_pass(all_cec_modules_df): + """ + Solve for all the modules in the DataFrame, using the empirical initial guess + """ + solver = pyo.SolverFactory('ipopt') + for i, r in all_cec_modules_df.iterrows(): if r['V_mp_ref'] < 0: all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" @@ -712,37 +567,57 @@ def solve_all_without_guess(all_cec_modules_df, plotting): all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" continue - model = create_model() + model = create_model(gamma_curve_dt=3) if not set_parameters(model.solver, r): - all_cec_modules_df.loc[i, 'Error'] = "Parameter missing or out of bounds" + all_cec_modules_df.loc[i, 'Error'] = f"Parameter missing or out of bounds for row {i}" continue if plotting: plt.figure() solver.options["max_iter"] = 3000 + set_empirical_initial_guess(model) + res, scaled_model = solve_model(model, solver, tee=False) - print(get_constraint_infeas(model)) if plotting: if res and res.solver.status == 'ok': plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) else: plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) - if res is None or res.solver.status != 'ok': + d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + + infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() + infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + + if infeas_max > 0.5 or infeas_sum > 0.5: + all_cec_modules_df.loc[i, 'Error'] = "Infeasibility > 0.5" + if plotting: + plt.close() continue a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) - all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] - all_cec_modules_df.to_csv("cec_modules_params.csv", index=False) + all_cec_modules_df.loc[i, iv_diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] + all_cec_modules_df.loc[i, model_param_cols] = [a, Il, Io, Rs, Rsh, Adj] + + if plotting: + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(plot_output_path / f"IV_curve_{i}.png") + plt.close() + return all_cec_modules_df -def solve_all_with_guess(all_cec_modules_df, plotting): - diff_cols = ['d_Isc', 'd_Imp', 'd_Vmp', 'd_Pmp'] +def run_solve_bootstrapping(all_cec_modules_df, solved_df, unsolved_df): + """ + Solve for all the modules in the DataFrame, using the `find_closest` to provide the initial guess + """ + solver = pyo.SolverFactory('ipopt') - unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] - solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] for i, r in unsolved_df.iterrows(): if r['V_mp_ref'] < 0: all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" @@ -751,9 +626,9 @@ def solve_all_with_guess(all_cec_modules_df, plotting): all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" continue - model = create_model(gamma_curve_dt=15) + model = create_model(gamma_curve_dt=3) if not set_parameters(model.solver, r): - all_cec_modules_df.loc[i, 'Error'] = "Parameter missing or out of bounds" + all_cec_modules_df.loc[i, 'Error'] = f"Parameter missing or out of bounds for row {i}" continue if plotting: @@ -761,13 +636,12 @@ def solve_all_with_guess(all_cec_modules_df, plotting): solver.options["max_iter"] = 3000 cec_closest = find_closest(solved_df, r) - set_initial_guess(model, *cec_closest[param_cols]) + set_initial_guess(model, *cec_closest[model_param_cols]) if plotting: plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") - res, scaled_model = solve_model_best_solution(model, solver, tee=False) + res, scaled_model = solve_model(model, solver, tee=False) - print(get_constraint_infeas(model)) if plotting: if res and res.solver.status == 'ok': plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) @@ -786,8 +660,9 @@ def solve_all_with_guess(all_cec_modules_df, plotting): continue a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) - all_cec_modules_df.loc[i, diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] - all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] + all_cec_modules_df.loc[i, iv_diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] + all_cec_modules_df.loc[i, model_param_cols] = [a, Il, Io, Rs, Rsh, Adj] + all_cec_modules_df.loc[i, 'Error'] = None if plotting: plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') @@ -795,32 +670,18 @@ def solve_all_with_guess(all_cec_modules_df, plotting): plt.ylabel("Current") plt.title(f"{r['Manufacturer']} {r['Model Number']}") plt.tight_layout() - plt.savefig(output_path / f"IV_curve_{'nosam' if not len(r_sam) else '_'}{i}.png") + plt.savefig(plot_output_path / f"IV_curve_{i}.png") plt.close() + return all_cec_modules_df - all_cec_modules_df.to_csv("cec_modules_params_approx.csv", index=False) - -if __name__ == "__main__": - - output_path = Path(__file__).parent / "6parsolve_output" - if not (output_path).exists(): - os.mkdir(output_path) +def run_solve_bootstrapping_reduced(all_cec_modules_df, solved_df, unsolved_df): + """ + Solve for all the modules in the DataFrame, using the closest initial guess and reducing the number of temperature samples to fit gamma + """ solver = pyo.SolverFactory('ipopt') - filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" - sam_solved_modules_df = pd.read_csv(filename, skiprows=[1, 2]) - - all_cec_modules_df = pd.read_csv("cec_modules_params_approx.csv") - - unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] - - solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] - - diff_cols = ['d_Isc', 'd_Imp', 'd_Vmp', 'd_Pmp'] - - plotting = True - for i, r in all_cec_modules_df.iterrows(): + for i, r in unsolved_df.iterrows(): if r['V_mp_ref'] < 0: all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" continue @@ -830,55 +691,26 @@ def solve_all_with_guess(all_cec_modules_df, plotting): model = create_model(gamma_curve_dt=15) if not set_parameters(model.solver, r): - all_cec_modules_df.loc[i, 'Error'] = "Parameter missing or out of bounds" + all_cec_modules_df.loc[i, 'Error'] = f"Parameter missing or out of bounds for row {i}" continue if plotting: plt.figure() solver.options["max_iter"] = 3000 - # plot SAM one - r_sam = sam_solved_modules_df[(sam_solved_modules_df['Manufacturer'] == r['Manufacturer'].replace(",", "")) - & (sam_solved_modules_df['Name'].str.contains(r['Model Number'], regex=False)) - & (sam_solved_modules_df['PTC'] == r['PTC'])].drop_duplicates() - if len(r_sam): - r_sam = r_sam.to_dict('records')[0] - set_initial_guess(model, r_sam['a_ref'], r_sam['I_L_ref'], r_sam['I_o_ref'], r_sam['R_s'], r_sam['R_sh_ref'], r_sam['Adjust']) - if plotting: - plot_iv_curve(model, linestyle='dashed', label="SAM") - if True: - cec_closest = find_closest(solved_df, r) - set_initial_guess(model, *cec_closest[param_cols]) - if plotting: - plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") + cec_closest = find_closest(solved_df, r) + set_initial_guess(model, *cec_closest[model_param_cols]) + if plotting: + plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") - # res, scaled_model = solve_model_best_solution(model, solver, tee=False) res, scaled_model = solve_model(model, solver, tee=False) - print(get_constraint_infeas(model)) if plotting: if res and res.solver.status == 'ok': plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) else: plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) - # if res is None or res.solver.status != 'ok': - # solver.options['max_iter'] = 1938 - # res = solver.solve(scaled_model, tee=True) - - # model_rel = create_full_relaxed_model() - # set_parameters(model_rel.solver, r) - - # copy_over_vars(model, scaled_model, model_rel) - - # res_rel, scaled_model_rel = solve_relaxed_model(model_rel, solver, True) - - # print(get_constraint_infeas(model_rel)) - # plot_iv_curve(model_rel, linestyle='-', label="Closest", plot_anchors=True) - - # if res_rel is None or res_rel.solver.status != 'ok': - # continue - d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() @@ -891,8 +723,9 @@ def solve_all_with_guess(all_cec_modules_df, plotting): continue a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) - all_cec_modules_df.loc[i, diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] - all_cec_modules_df.loc[i, param_cols] = [a, Il, Io, Rs, Rsh, Adj] + all_cec_modules_df.loc[i, iv_diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] + all_cec_modules_df.loc[i, model_param_cols] = [a, Il, Io, Rs, Rsh, Adj] + all_cec_modules_df.loc[i, 'Error'] = None if plotting: plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') @@ -900,373 +733,142 @@ def solve_all_with_guess(all_cec_modules_df, plotting): plt.ylabel("Current") plt.title(f"{r['Manufacturer']} {r['Model Number']}") plt.tight_layout() - plt.savefig(output_path / f"IV_curve_{'nosam' if not len(r_sam) else '_'}{i}.png") + plt.savefig(plot_output_path / f"IV_curve_{i}.png") plt.close() + return all_cec_modules_df - all_cec_modules_df.to_csv("cec_modules_params_approx.csv", index=False) - exit() - - - - - - - filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" - cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) - - filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/Bad Modules 2024-Nov-15_09-53.csv" - cec_bad_modules_df = pd.read_csv(filename, skiprows=[1, 2]) - cec_bad_modules_df.index = range(len(cec_modules_df), len(cec_modules_df) + len(cec_bad_modules_df)) - - cec_bad_modules_df = cec_bad_modules_df.rename(columns={"Model": "Name", "celltype": "Technology", "vmp": "V_mp_ref", - "imp": "I_mp_ref", "voc": "V_oc_ref", "isc": "I_sc_ref", "alpha_isc": "alpha_sc", - "beta_voc": "beta_oc", "gamma_pmp": "gamma_r", "n_ser": "N_s", "t_ref": "STC"}) - cec_bad_modules_df = cec_bad_modules_df.drop(columns=["Message"]) - - - solved_modules = "/Users/dguittet/Projects/SAM/pysam/param_comparison_1.csv" - solved_modules_df = pd.read_csv(solved_modules, index_col='index').dropna() - cec_modules_solved_df = cec_modules_df[cec_modules_df.index.isin(solved_modules_df.index)] - solved_bad_modules = "/Users/dguittet/Projects/SAM/pysam/bad_modules_param_0.csv" - solved_bad_modules_df = pd.read_csv(solved_bad_modules, index_col='index') - solved_bad_modules_df.index = range(len(cec_modules_df), len(cec_modules_df) + len(solved_bad_modules_df)) - cec_bad_modules_solved_df = cec_bad_modules_df[cec_bad_modules_df.index.isin(solved_bad_modules_df.index)] +def solve_bootstrapping_reduced_approx(all_cec_modules_df, solved_df, unsolved_df): + """ + Solve for all the modules in the DataFrame, using the closest initial guess, reducing the number of temperature samples to fit gamma, + and using `solve_model_best_solution` - # read all - filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/PV_Module_List_Full_Data_ADA-2024-11-12.xlsx" - all_cec_modules_df = pd.read_excel(filename, skiprows=list(range(0, 16)) + [17]) - all_cec_modules_df = all_cec_modules_df.drop(columns=["Description", 'Safety Certification', - 'Nameplate Pmax', 'Notes', - 'Design Qualification Certification\n(Optional Submission)', - 'Performance Evaluation (Optional Submission)', 'Family', - 'N_p', 'αIpmax', 'βVpmax', 'IPmax, low', 'VPmax, low', 'IPmax, NOCT', - 'VPmax, NOCT', 'Mounting', 'Type', 'Short Side', 'Long Side', - 'Geometric Multiplier', 'P2/Pref', 'CEC Listing Date', 'Last Update']) - all_cec_modules_df = all_cec_modules_df.rename(columns={ - 'Nameplate Isc': "I_sc_ref", 'Nameplate Voc': "V_oc_ref", - 'Nameplate Ipmax': "I_mp_ref", 'Nameplate Vpmax': "V_mp_ref", 'Average NOCT': "T_NOCT", - 'γPmax': "gamma_r", 'αIsc': "alpha_sc", - 'βVoc': "beta_oc", - }) - num_cols = ['A_c', 'N_s', 'I_sc_ref', 'V_oc_ref', 'I_mp_ref', 'V_mp_ref', 'T_NOCT', - 'gamma_r', 'alpha_sc', 'beta_oc'] - all_cec_modules_df[num_cols] = all_cec_modules_df[num_cols].astype(float) - all_cec_modules_df['alpha_sc'] *= all_cec_modules_df['I_sc_ref'] * 1e-2 - all_cec_modules_df['beta_oc'] *= all_cec_modules_df['V_oc_ref'] * 1e-2 - all_cec_modules_df = all_cec_modules_df.drop_duplicates() - - for col in param_cols: - all_cec_modules_df[col] = None + Not recommended to use this function unless an approximate solution is desperately needed. + Accuracy of IV curve to test data should be examined visually. This can be done with `create_model_with_solution` and `plot_iv_curve` + """ + solver = pyo.SolverFactory('ipopt') - for i, row in all_cec_modules_df.iterrows(): - r_sam = cec_modules_solved_df[(cec_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) - & (cec_modules_solved_df['Name'].str.contains(row['Model Number'], regex=False)) - & (cec_modules_solved_df['PTC'] == row['PTC'])].drop_duplicates() - if len(r_sam) > 1: - r_sam = r_sam[r_sam['Name'].str.endswith(row['Model Number'])] - if len(r_sam) > 1: - print("Duplicate: ", i, row) - continue - if len(r_sam) == 1: - param_row = solved_modules_df.loc[r_sam.index] - - if len(r_sam) == 0: - r_sam = cec_bad_modules_solved_df[(cec_bad_modules_solved_df['Manufacturer'] == row['Manufacturer'].replace(",", "")) - & (cec_bad_modules_solved_df['Name'].str.contains(row['Model Number'], regex=False)) - ].drop_duplicates() - if len(r_sam) > 1: - r_sam = r_sam[r_sam['Name'].str.endswith(row['Model Number'])] - if len(r_sam) > 1: - print("Duplicate: ", i, row) - continue - if len(r_sam) == 1: - param_row = solved_bad_modules_df.loc[r_sam.index] - - if len(r_sam) == 0: + for i, r in unsolved_df.iterrows(): + if r['V_mp_ref'] < 0: + all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" + continue + if r['V_oc_ref'] < r['V_mp_ref']: + all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" continue - all_cec_modules_df.loc[i, param_cols] = param_row[param_cols].values[0] - all_cec_modules_df.to_csv("cec_modules_params.csv") - - exit() - + model = create_model(gamma_curve_dt=15) + if not set_parameters(model.solver, r): + all_cec_modules_df.loc[i, 'Error'] = f"Parameter missing or out of bounds for row {i}" + continue + if plotting: + plt.figure() + solver.options["max_iter"] = 3000 + cec_closest = find_closest(solved_df, r) + set_initial_guess(model, *cec_closest[model_param_cols]) + if plotting: + plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") + res, scaled_model = solve_model_best_solution(model, solver, tee=False) + + if plotting: + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) + else: + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) - param_comparison = { - "index": [], - "a_ssc": [], - "Il_ssc": [], - "Io_ssc": [], - "Rs_ssc": [], - "Rsh_ssc": [], - "Adj_ssc": [], - "a_py": [], - "Il_py": [], - "Io_py": [], - "Rs_py": [], - "Rsh_py": [], - "Adj_py": [] - } + d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) - for i, r in cec_modules_unsolved_df.iterrows(): - model = create_model() - m = model.solver - tech = r['Technology'] + infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() + infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + if infeas_max > 0.5 or infeas_sum > 0.5: + all_cec_modules_df.loc[i, 'Error'] = "Infeasibility > 0.5" + if plotting: + plt.close() + continue - param_comparison['index'].append(i) - param_comparison['a_ssc'].append(r['a_ref']) - param_comparison['Il_ssc'].append(r['I_L_ref']) - param_comparison['Io_ssc'].append(r['I_o_ref']) - param_comparison['Rs_ssc'].append(r['R_s']) - param_comparison['Rsh_ssc'].append(r['R_sh_ref']) - param_comparison['Adj_ssc'].append(r['Adjust']) + a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) + all_cec_modules_df.loc[i, iv_diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] + all_cec_modules_df.loc[i, model_param_cols] = [a, Il, Io, Rs, Rsh, Adj] + all_cec_modules_df.loc[i, 'Error'] = None - m.Vmp.set_value(r['V_mp_ref']) - m.Imp.set_value(r['I_mp_ref']) - m.Voc.set_value(r['V_oc_ref']) - m.Isc.set_value(r['I_sc_ref']) - m.aIsc.set_value(r['alpha_sc']) - m.bVoc.set_value(r['beta_oc']) - m.gPmp.set_value(r['gamma_r']) - # m.Tref.set_value(r['STC']) - m.Tref.set_value(25 + 273.15) + if plotting: + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(plot_output_path / f"IV_curve_{i}.png") + plt.close() + return all_cec_modules_df - df = cec_modules_solved_df - diff = ( - (df['V_mp_ref'] - r['V_mp_ref']) ** 2 + - (df['I_mp_ref'] - r['I_mp_ref']) ** 2 + - (df['V_oc_ref'] - r['V_oc_ref']) ** 2 + - (df['I_sc_ref'] - r['I_sc_ref']) ** 2 + - (df['alpha_sc'] - r['alpha_sc']) ** 2 + - (df['beta_oc'] - r['beta_oc']) ** 2 + - (df['gamma_r'] - r['gamma_r']) ** 2) - - cec_closest = cec_modules_solved_df.iloc[diff.argmin()] - m.par.a.set_value(cec_closest['a_ref']) - m.par.Il.set_value(cec_closest['I_L_ref']) - m.par.Io.set_value(cec_closest['I_o_ref']) - m.par.Rs.set_value(cec_closest['R_s']) - m.par.Rsh.set_value(cec_closest['R_sh_ref']) - m.par.Adj.set_value(cec_closest['Adjust']) +def create_model_with_solution(sol_row): + """ + Create a model with the parameters from the tests and for the single-diode model. Can then be used for `plot_iv_curve` + + sol_row must contain columns from `test_data_cols` and `model_param_cols` + """ + model = create_model(gamma_curve_dt=3) + if type(sol_row) == pd.Series: + if not set_parameters(model.solver, sol_row): + raise RuntimeError("Test data parameters could not be set") + set_initial_guess(model, *sol_row[model_param_cols]) + elif type(sol_row) == pd.DataFrame: + if not set_parameters(model.solver, sol_row.to_dict('records')[0]): + raise RuntimeError("Test data parameters could not be set") + set_initial_guess(model, *sol_row[model_param_cols].values[0]) + return model - res, scaled_model = solve_model(model, solver, tee=False) - if res is None or res.solver.status != 'ok': - # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) - # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) - # print(pyo.value(scaled_model.solver.gamma.obj)) - param_comparison['a_py'].append(None) - param_comparison['Il_py'].append(None) - param_comparison['Io_py'].append(None) - param_comparison['Rs_py'].append(None) - param_comparison['Rsh_py'].append(None) - param_comparison['Adj_py'].append(None) - else: +if __name__ == "__main__": + filename = None - a = pyo.value(scaled_model.solver.par.scaled_a) - Il = pyo.value(scaled_model.solver.par.scaled_Il) - Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) - Rs = pyo.value(scaled_model.solver.par.scaled_Rs) - Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) - Adj = pyo.value(scaled_model.solver.par.scaled_Adj) - - param_comparison['a_py'].append(a) - param_comparison['Il_py'].append(Il) - param_comparison['Io_py'].append(Io) - param_comparison['Rs_py'].append(Rs) - param_comparison['Rsh_py'].append(Rsh) - param_comparison['Adj_py'].append(Adj) - - # if i % 100 == 0: - # comparison_df = pd.DataFrame(param_comparison).set_index('index') - # comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() - # comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() - # comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() - # comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() - # comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() - # comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() - - # comparison_df = comparison_df.dropna() - # cec_modules_solved_df = pd.concat([cec_modules_solved_df, comparison_df]) - - comparison_df = pd.DataFrame(param_comparison).set_index('index') - comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() - comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() - comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() - comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() - comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() - comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() - - comparison_df = comparison_df.dropna() - for i, row in comparison_df.iterrows(): - solved_modules_df.loc[i] = row - solved_modules_df.to_csv("param_comparison.csv") - - exit() - filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/Bad Modules 2024-Nov-15_09-53.csv" - bad_modules_df = pd.read_csv(filename, skiprows=[1, 2]) - - param_comparison = { - "index": [], - "a_py": [], - "Il_py": [], - "Io_py": [], - "Rs_py": [], - "Rsh_py": [], - "Adj_py": [] - } - - tee = False - - for n, r in bad_modules_df.iterrows(): - model = create_model() - m = model.solver - tech = r['Technology'] - - param_comparison['index'].append(n) - - m.Vmp.set_value(r['vmp']) - m.Imp.set_value(r['imp']) - m.Voc.set_value(r['voc']) - m.Isc.set_value(r['isc']) - m.aIsc.set_value(r['alpha_isc']) - m.bVoc.set_value(r['beta_voc']) - m.gPmp.set_value(r['gamma_pmp']) - m.Tref.set_value(25 + 273.15) + if len(sys.argv) > 1: + filename = Path(sys.argv[1]) - res, scaled_model = solve_model(model, solver, tee=tee) - - if res is None or res.solver.status != 'ok': - # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) - # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) - # print(pyo.value(scaled_model.solver.gamma.obj)) - param_comparison['a_py'].append(None) - param_comparison['Il_py'].append(None) - param_comparison['Io_py'].append(None) - param_comparison['Rs_py'].append(None) - param_comparison['Rsh_py'].append(None) - param_comparison['Adj_py'].append(None) - else: + if not filename.exists(): + raise RuntimeError(f"CEC Module Excel Spreadsheet file path does not exist. {filename}") + + filename_date = filename.stem.split("-") + filename_date = f"{filename_date[-3]}-{filename_date[-2]}-{filename_date[-1]}" + all_cec_modules_df = read_prepare_file(filename) - a = pyo.value(scaled_model.solver.par.scaled_a) - Il = pyo.value(scaled_model.solver.par.scaled_Il) - Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) - Rs = pyo.value(scaled_model.solver.par.scaled_Rs) - Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) - Adj = pyo.value(scaled_model.solver.par.scaled_Adj) + plotting = False - param_comparison['a_py'].append(a) - param_comparison['Il_py'].append(Il) - param_comparison['Io_py'].append(Io) - param_comparison['Rs_py'].append(Rs) - param_comparison['Rsh_py'].append(Rsh) - param_comparison['Adj_py'].append(Adj) + if plotting: + plot_output_path = Path(__file__).parent / "6parsolve_output" + if not (plot_output_path).exists(): + os.mkdir(plot_output_path) - if n % 10 == 0: - comparison_df = pd.DataFrame(param_comparison).set_index('index') - comparison_df.to_csv("bad_modules_param.csv") + print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Solving for {len(all_cec_modules_df)} Modules") + # solve attempt #1 + print("Starting First Pass Solve") + all_cec_modules_df = run_solve_first_pass(all_cec_modules_df) + all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) + + solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] + unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] + print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Solved {len(solved_df)}. {len(unsolved_df)} remaining.") -def run_solved_modules(): - solver = pyo.SolverFactory('ipopt') - solver.options["max_iter"] = 5000 - - filename = "/Users/dguittet/Projects/SAM/sam/samples/CEC Module and Inverter Libraries/CEC Modules/CEC Modules 2024-11-14/CEC Modules 2024-Nov-15_09-53.csv" - cec_modules_df = pd.read_csv(filename, skiprows=[1, 2]) - - param_comparison = { - "index": [], - "a_ssc": [], - "Il_ssc": [], - "Io_ssc": [], - "Rs_ssc": [], - "Rsh_ssc": [], - "Adj_ssc": [], - "a_py": [], - "Il_py": [], - "Io_py": [], - "Rs_py": [], - "Rsh_py": [], - "Adj_py": [] - } - - tee = False - - for n, r in cec_modules_df.iterrows(): - model = create_model() - m = model.solver - tech = r['Technology'] - - # if r['I_mp_ref'] < r['I_sc_ref']: - # raise ValueError - - param_comparison['index'].append(n) - param_comparison['a_ssc'].append(r['a_ref']) - param_comparison['Il_ssc'].append(r['I_L_ref']) - param_comparison['Io_ssc'].append(r['I_o_ref']) - param_comparison['Rs_ssc'].append(r['R_s']) - param_comparison['Rsh_ssc'].append(r['R_sh_ref']) - param_comparison['Adj_ssc'].append(r['Adjust']) + # solve with bootstrapping + print("{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Second Pass Solve") + all_cec_modules_df = run_solve_bootstrapping(all_cec_modules_df, solved_df, unsolved_df) + all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) - m.Vmp.set_value(r['V_mp_ref']) - m.Imp.set_value(r['I_mp_ref']) - m.Voc.set_value(r['V_oc_ref']) - m.Isc.set_value(r['I_sc_ref']) - m.aIsc.set_value(r['alpha_sc']) - m.bVoc.set_value(r['beta_oc']) - m.gPmp.set_value(r['gamma_r']) - m.Tref.set_value(25 + 273.15) + solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] + unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] + print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Solved {len(solved_df)}. {len(unsolved_df)} remaining.") - # m.par.a.set_value(r['a_ref']) - # m.par.Il.set_value(r['I_L_ref']) - # m.par.Io.set_value(r['I_o_ref']) - # m.par.Rs.set_value(r['R_s']) - # m.par.Rsh.set_value(r['R_sh_ref']) - # m.par.Adj.set_value(r['Adjust']) - - res, scaled_model = solve_model(model, solver, tee=tee) - - if res is None or res.solver.status != 'ok': - # log_infeasible_constraints(scaled_model, logger=solve_log, tol=1e-5, log_expression=True, log_variables=True) - # log_infeasible_bounds(scaled_model, logger=solve_log, tol=1e-5) - # print(n, pyo.value(scaled_model.solver.gamma.gamma_avg), pyo.value(scaled_model.solver.gPmp)) - # print(pyo.value(scaled_model.solver.gamma.obj)) - param_comparison['a_py'].append(None) - param_comparison['Il_py'].append(None) - param_comparison['Io_py'].append(None) - param_comparison['Rs_py'].append(None) - param_comparison['Rsh_py'].append(None) - param_comparison['Adj_py'].append(None) - else: + # solve with bootstrapping & reduced number of temperature samples + print("{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Final Pass Solve") + all_cec_modules_df = run_solve_bootstrapping_reduced(all_cec_modules_df, solved_df, unsolved_df) + all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) - a = pyo.value(scaled_model.solver.par.scaled_a) - Il = pyo.value(scaled_model.solver.par.scaled_Il) - Io = pyo.value(scaled_model.solver.par.scaled_Io / IL_SCALING) - Rs = pyo.value(scaled_model.solver.par.scaled_Rs) - Rsh = pyo.value(scaled_model.solver.par.scaled_Rsh / RSH_SCALING) - Adj = pyo.value(scaled_model.solver.par.scaled_Adj) - - param_comparison['a_py'].append(a) - param_comparison['Il_py'].append(Il) - param_comparison['Io_py'].append(Io) - param_comparison['Rs_py'].append(Rs) - param_comparison['Rsh_py'].append(Rsh) - param_comparison['Adj_py'].append(Adj) - - if n % 10 == 0: - comparison_df = pd.DataFrame(param_comparison).set_index('index') - comparison_df['d_a'] = ((comparison_df['a_py'] - comparison_df['a_ssc'])/comparison_df['a_ssc']).abs() - comparison_df['d_Il'] = ((comparison_df['Il_py'] - comparison_df['Il_ssc'])/comparison_df['Il_ssc']).abs() - comparison_df['d_Io'] = ((comparison_df['Io_py'] - comparison_df['Io_ssc'])/comparison_df['Io_ssc']).abs() - comparison_df['d_Rs'] = ((comparison_df['Rs_py'] - comparison_df['Rs_ssc'])/comparison_df['Rs_ssc']).abs() - comparison_df['d_Rsh'] = ((comparison_df['Rsh_py'] - comparison_df['Rsh_ssc'])/comparison_df['Rsh_ssc']).abs() - comparison_df['d_Adj'] = ((comparison_df['Adj_py'] - comparison_df['Adj_ssc'])/comparison_df['Adj_ssc']).abs() - - comparison_df.to_csv("param_comparison.csv") + solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] + unsolved_df = all_cec_modules_df[all_cec_modules_df['Rsh_py'].isna()] + print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Solved {len(solved_df)}. {len(unsolved_df)} unsolved.") + \ No newline at end of file From 54dadca4420a2711d0ab39c273a618195ee6d4dd Mon Sep 17 00:00:00 2001 From: Darice Date: Tue, 27 May 2025 12:42:38 -0600 Subject: [PATCH 09/17] update for tests --- .github/workflows/test_pkg.yml | 3 +++ tests/requirements.txt | 1 + tests/test_SixParSolve.py | 6 ------ 3 files changed, 4 insertions(+), 6 deletions(-) diff --git a/.github/workflows/test_pkg.yml b/.github/workflows/test_pkg.yml index 82b0179..539f5b1 100644 --- a/.github/workflows/test_pkg.yml +++ b/.github/workflows/test_pkg.yml @@ -38,6 +38,7 @@ jobs: pip install -r requirements.txt pip install -r tests/requirements.txt pip install nrel-pysam==${{ env.VERSION }} + conda install -c conda-forge ipopt - uses: actions/checkout@v3 - name: Checkout SAM @@ -89,6 +90,7 @@ jobs: pip install -r requirements.txt pip install -r tests/requirements.txt pip install nrel-pysam==$VER + conda install -c conda-forge ipopt - uses: actions/checkout@v3 - name: Checkout SAM @@ -142,6 +144,7 @@ jobs: pip install -r requirements.txt pip install -r tests/requirements.txt pip install nrel-pysam==$VER + conda install -c conda-forge ipopt - uses: actions/checkout@v3 - name: Checkout SAM diff --git a/tests/requirements.txt b/tests/requirements.txt index 038d4ac..6dd8c17 100644 --- a/tests/requirements.txt +++ b/tests/requirements.txt @@ -1,6 +1,7 @@ mypy numpy pandas +pyomo pytest python-dotenv pympler diff --git a/tests/test_SixParSolve.py b/tests/test_SixParSolve.py index cbc15cd..f98925c 100644 --- a/tests/test_SixParSolve.py +++ b/tests/test_SixParSolve.py @@ -96,9 +96,3 @@ def test_cec_model_solve(set_initial_values=True): Adj = pyo.value(scaled_model.solver.par.scaled_Adj) plot_iv_curve(model) - -# fig = plt.figure() -# test_default_sam_cec_user() - -# test_cec_model_solve(set_initial_values=False) -# plt.show() \ No newline at end of file From 926a311294c7c72060294048824895cd09e434dd Mon Sep 17 00:00:00 2001 From: Darice Date: Tue, 27 May 2025 12:45:16 -0600 Subject: [PATCH 10/17] update test requirements.txt --- tests/requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/requirements.txt b/tests/requirements.txt index 6dd8c17..a179afa 100644 --- a/tests/requirements.txt +++ b/tests/requirements.txt @@ -1,3 +1,4 @@ +matplotlib mypy numpy pandas From 18228b29050ee21ce4679f3d5e91d3fb34b3e3cb Mon Sep 17 00:00:00 2001 From: Darice Date: Tue, 27 May 2025 12:57:32 -0600 Subject: [PATCH 11/17] fix tests --- .github/workflows/test_pkg.yml | 1 - files/SixParSolve.py | 5 ++++- tests/test_SixParSolve.py | 23 ++++++++++++----------- 3 files changed, 16 insertions(+), 13 deletions(-) diff --git a/.github/workflows/test_pkg.yml b/.github/workflows/test_pkg.yml index 539f5b1..ff8fc60 100644 --- a/.github/workflows/test_pkg.yml +++ b/.github/workflows/test_pkg.yml @@ -38,7 +38,6 @@ jobs: pip install -r requirements.txt pip install -r tests/requirements.txt pip install nrel-pysam==${{ env.VERSION }} - conda install -c conda-forge ipopt - uses: actions/checkout@v3 - name: Checkout SAM diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 3af4d26..572e4fd 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -11,7 +11,10 @@ from datetime import datetime logging.getLogger('pyomo.core').setLevel(logging.ERROR) -solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties") +solve_log = logging.getLogger('solve_log') +solve_log.setLevel(logging.INFO) +solve_log = logging.LoggerAdapter(solve_log, {"tag": solve_log}) + IL_SCALING = 1e8 RSH_SCALING = 1e-3 test_data_cols = ['A_c', 'N_s', 'I_sc_ref', 'V_oc_ref', 'I_mp_ref', 'V_mp_ref', 'T_NOCT', diff --git a/tests/test_SixParSolve.py b/tests/test_SixParSolve.py index f98925c..07083de 100644 --- a/tests/test_SixParSolve.py +++ b/tests/test_SixParSolve.py @@ -1,4 +1,4 @@ -from pytest import approx +import pytest import sys from pathlib import Path sys.path.append(str(Path(__file__).parent.parent)) @@ -46,11 +46,11 @@ def test_default_sam_cec_user(): model = define_default_model() IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(model, 1000, 25+275.15) - assert IL_oper == approx(6.059759, rel=1e-5) - assert IO_oper == approx(1.1674203993060455e-10, rel=1e-5) - assert Rs == approx(0.3081202, rel=1e-5) - assert A_oper == approx(2.5948906, rel=1e-5) - assert Rsh_oper == approx(500.069, rel=1e-5) + assert IL_oper == pytest.approx(6.059759, rel=1e-5) + assert IO_oper == pytest.approx(1.1674203993060455e-10, rel=1e-5) + assert Rs == pytest.approx(0.3081202, rel=1e-5) + assert A_oper == pytest.approx(2.5948906, rel=1e-5) + assert Rsh_oper == pytest.approx(500.069, rel=1e-5) plot_iv_curve(model) @@ -73,13 +73,14 @@ def test_cec_model_ivcurve_default(): I = current_at_voltage_cec( v, IL_oper, IO_oper, Rs, A_oper, Rsh_oper, I_mp_ref ) y_I.append(I) - assert V[1] == approx(0.43221, rel=1e-3) - assert V[-1] == approx(64.4, rel=1e-3) - assert y_I[0] == approx(6.05, rel=1e-3) - assert y_I[1] == approx(6.04913, rel=1e-3) - assert y_I[-1] == approx(5.770225265737295e-06, rel=1e-3) + assert V[1] == pytest.approx(0.43221, rel=1e-3) + assert V[-1] == pytest.approx(64.4, rel=1e-3) + assert y_I[0] == pytest.approx(6.05, rel=1e-3) + assert y_I[1] == pytest.approx(6.04913, rel=1e-3) + assert y_I[-1] == pytest.approx(5.770225265737295e-06, rel=1e-3) +@pytest.mark.skipif(not pyo.SolverFactory('ipopt').available(), reason="requires ipopt solver") def test_cec_model_solve(set_initial_values=True): model = define_default_model(set_initial_values) From 85a15da972748b82399fc46e9201e6fd1fd28ec9 Mon Sep 17 00:00:00 2001 From: Darice Date: Tue, 27 May 2025 13:01:03 -0600 Subject: [PATCH 12/17] fix imports --- files/SixParSolve.py | 1 - 1 file changed, 1 deletion(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 572e4fd..8e2829d 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -5,7 +5,6 @@ import pandas as pd import numpy as np import matplotlib.pyplot as plt -import idaes.logger as idaeslog from pyomo.util.infeasible import log_infeasible_constraints, log_infeasible_bounds, log_close_to_bounds import logging from datetime import datetime From 273c8d95c78ea9c309b269fd9caae7c73df2dca3 Mon Sep 17 00:00:00 2001 From: Darice Date: Tue, 27 May 2025 13:06:45 -0600 Subject: [PATCH 13/17] update requirements.txt --- tests/requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/requirements.txt b/tests/requirements.txt index a179afa..f6c9c60 100644 --- a/tests/requirements.txt +++ b/tests/requirements.txt @@ -2,6 +2,7 @@ matplotlib mypy numpy pandas +pint pyomo pytest python-dotenv From 56c1d30dba2909c20e30753e7b4b5d98fb7e159a Mon Sep 17 00:00:00 2001 From: Darice Date: Wed, 28 May 2025 08:04:48 -0600 Subject: [PATCH 14/17] add multiprocessing --- files/SixParSolve.py | 352 +++++++++++++++++++++++++------------------ 1 file changed, 207 insertions(+), 145 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 8e2829d..36028d9 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -8,6 +8,9 @@ from pyomo.util.infeasible import log_infeasible_constraints, log_infeasible_bounds, log_close_to_bounds import logging from datetime import datetime +import multiprocessing as mp + +NUM_OF_WORKERS = mp.cpu_count() // 2 logging.getLogger('pyomo.core').setLevel(logging.ERROR) solve_log = logging.getLogger('solve_log') @@ -555,189 +558,246 @@ def read_prepare_file(xlsx_file): return all_cec_modules_df -def run_solve_first_pass(all_cec_modules_df): +def run_solve_first_pass(i, r, plotting=False, solver=None): """ Solve for all the modules in the DataFrame, using the empirical initial guess """ - solver = pyo.SolverFactory('ipopt') + if solver is None: + solver = pyo.SolverFactory('ipopt') - for i, r in all_cec_modules_df.iterrows(): - if r['V_mp_ref'] < 0: - all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" - continue - if r['V_oc_ref'] < r['V_mp_ref']: - all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" - continue + if r['V_mp_ref'] < 0: + r['Error'] = "Vmp < 0" + return r + if r['V_oc_ref'] < r['V_mp_ref']: + r['Error'] = "Voc < Vmp" + return r - model = create_model(gamma_curve_dt=3) - if not set_parameters(model.solver, r): - all_cec_modules_df.loc[i, 'Error'] = f"Parameter missing or out of bounds for row {i}" - continue - - if plotting: - plt.figure() - solver.options["max_iter"] = 3000 - - set_empirical_initial_guess(model) + model = create_model(gamma_curve_dt=3) + if not set_parameters(model.solver, r): + r['Error'] = f"Parameter missing or out of bounds for row {i}" + return r - res, scaled_model = solve_model(model, solver, tee=False) - - if plotting: - if res and res.solver.status == 'ok': - plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) - else: - plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) + if plotting: + plt.figure() + solver.options["max_iter"] = 3000 - d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + set_empirical_initial_guess(model) - infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() - infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + res, scaled_model = solve_model(model, solver, tee=False) + + if plotting: + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) + else: + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) - if infeas_max > 0.5 or infeas_sum > 0.5: - all_cec_modules_df.loc[i, 'Error'] = "Infeasibility > 0.5" - if plotting: - plt.close() - continue + d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) - a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) - all_cec_modules_df.loc[i, iv_diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] - all_cec_modules_df.loc[i, model_param_cols] = [a, Il, Io, Rs, Rsh, Adj] + infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() + infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + if infeas_max > 0.5 or infeas_sum > 0.5: + r['Error'] = "Infeasibility > 0.5" if plotting: - plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - plt.xlabel("Voltage") - plt.ylabel("Current") - plt.title(f"{r['Manufacturer']} {r['Model Number']}") - plt.tight_layout() - plt.savefig(plot_output_path / f"IV_curve_{i}.png") plt.close() - return all_cec_modules_df + return r + a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) + for col, val in zip(iv_diff_cols, [d_I_sc, d_I_mp, d_V_mp, d_P_mp]): + r[col] = val + for col, val in zip(model_param_cols, [a, Il, Io, Rs, Rsh, Adj]): + r[col] = val -def run_solve_bootstrapping(all_cec_modules_df, solved_df, unsolved_df): - """ - Solve for all the modules in the DataFrame, using the `find_closest` to provide the initial guess - """ + if plotting: + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(plot_output_path / f"IV_curve_{i}.png") + plt.close() + return r + + +def parallel_run_solve_first_pass(all_cec_modules_df, plotting): + with mp.Pool(NUM_OF_WORKERS) as pool: + results = [pool.apply_async(run_solve_first_pass, [idx, row, plotting]) for idx, row in all_cec_modules_df.iterrows()] + results = [row.get() for row in results] + df = pd.DataFrame(results) + return df + +def sequential_run_solve_first_pass(all_cec_modules_df, plotting=False): solver = pyo.SolverFactory('ipopt') - for i, r in unsolved_df.iterrows(): - if r['V_mp_ref'] < 0: - all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" - continue - if r['V_oc_ref'] < r['V_mp_ref']: - all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" - continue + results = [] + for i, r in all_cec_modules_df.iterrows(): + row = run_solve_first_pass(i, r, plotting, solver) + results.append(row) + return pd.DataFrame(results) - model = create_model(gamma_curve_dt=3) - if not set_parameters(model.solver, r): - all_cec_modules_df.loc[i, 'Error'] = f"Parameter missing or out of bounds for row {i}" - continue - if plotting: - plt.figure() - solver.options["max_iter"] = 3000 +def run_solve_bootstrapping(i, r, solved_df, plotting=False, solver=None): + """ + Solve for all the modules in the DataFrame, using the `find_closest` to provide the initial guess + """ + if solver is None: + solver = pyo.SolverFactory('ipopt') + solver.options["max_iter"] = 3000 - cec_closest = find_closest(solved_df, r) - set_initial_guess(model, *cec_closest[model_param_cols]) - if plotting: - plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") + if r['V_mp_ref'] < 0: + r['Error'] = "Vmp < 0" + return r + if r['V_oc_ref'] < r['V_mp_ref']: + r['Error'] = "Voc < Vmp" + return r - res, scaled_model = solve_model(model, solver, tee=False) - - if plotting: - if res and res.solver.status == 'ok': - plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) - else: - plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) + model = create_model(gamma_curve_dt=3) + if not set_parameters(model.solver, r): + r['Error'] = f"Parameter missing or out of bounds for row {i}" + return r - d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + if plotting: + plt.figure() - infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() - infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + cec_closest = find_closest(solved_df, r) + set_initial_guess(model, *cec_closest[model_param_cols]) + if plotting: + plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") - if infeas_max > 0.5 or infeas_sum > 0.5: - all_cec_modules_df.loc[i, 'Error'] = "Infeasibility > 0.5" - if plotting: - plt.close() - continue + res, scaled_model = solve_model(model, solver, tee=False) + + if plotting: + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) + else: + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) - a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) - all_cec_modules_df.loc[i, iv_diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] - all_cec_modules_df.loc[i, model_param_cols] = [a, Il, Io, Rs, Rsh, Adj] - all_cec_modules_df.loc[i, 'Error'] = None + d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + + infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() + infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + if infeas_max > 0.5 or infeas_sum > 0.5: + r['Error'] = "Infeasibility > 0.5" if plotting: - plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - plt.xlabel("Voltage") - plt.ylabel("Current") - plt.title(f"{r['Manufacturer']} {r['Model Number']}") - plt.tight_layout() - plt.savefig(plot_output_path / f"IV_curve_{i}.png") plt.close() - return all_cec_modules_df + return r + a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) + for col, val in zip(iv_diff_cols, [d_I_sc, d_I_mp, d_V_mp, d_P_mp]): + r[col] = val + for col, val in zip(model_param_cols, [a, Il, Io, Rs, Rsh, Adj]): + r[col] = val + r['Error'] = None -def run_solve_bootstrapping_reduced(all_cec_modules_df, solved_df, unsolved_df): - """ - Solve for all the modules in the DataFrame, using the closest initial guess and reducing the number of temperature samples to fit gamma - """ + if plotting: + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(plot_output_path / f"IV_curve_{i}.png") + plt.close() + return r + + +def parallel_run_solve_bootstrapping(solved_df, unsolved_df, plotting=False): + with mp.Pool(NUM_OF_WORKERS) as pool: + results = [pool.apply_async(run_solve_bootstrapping, [idx, row, solved_df, plotting]) for idx, row in unsolved_df.iterrows()] + results = [row.get() for row in results] + df = pd.DataFrame(results) + return df + +def sequential_run_solve_bootstrapping(solved_df, unsolved_df, plotting=False): solver = pyo.SolverFactory('ipopt') + results = [] for i, r in unsolved_df.iterrows(): - if r['V_mp_ref'] < 0: - all_cec_modules_df.loc[i, 'Error'] = "Vmp < 0" - continue - if r['V_oc_ref'] < r['V_mp_ref']: - all_cec_modules_df.loc[i, 'Error'] = "Voc < Vmp" - continue + row = run_solve_bootstrapping(i, r, solved_df, plotting, solver) + results.append(row) + return pd.DataFrame(results) - model = create_model(gamma_curve_dt=15) - if not set_parameters(model.solver, r): - all_cec_modules_df.loc[i, 'Error'] = f"Parameter missing or out of bounds for row {i}" - continue - if plotting: - plt.figure() - solver.options["max_iter"] = 3000 +def run_solve_bootstrapping_reduced(i, r, solved_df, plotting=False, solver=None): + """ + Solve for all the modules in the DataFrame, using the closest initial guess and reducing the number of temperature samples to fit gamma + """ + if solver is None: + solver = pyo.SolverFactory('ipopt') + solver.options["max_iter"] = 3000 - cec_closest = find_closest(solved_df, r) - set_initial_guess(model, *cec_closest[model_param_cols]) - if plotting: - plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") + if r['V_mp_ref'] < 0: + r['Error'] = "Vmp < 0" + return r + if r['V_oc_ref'] < r['V_mp_ref']: + r['Error'] = "Voc < Vmp" + return r - res, scaled_model = solve_model(model, solver, tee=False) - - if plotting: - if res and res.solver.status == 'ok': - plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) - else: - plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) + model = create_model(gamma_curve_dt=15) + if not set_parameters(model.solver, r): + r['Error'] = f"Parameter missing or out of bounds for row {i}" + return r - d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + if plotting: + plt.figure() - infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() - infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + cec_closest = find_closest(solved_df, r) + set_initial_guess(model, *cec_closest[model_param_cols]) + if plotting: + plot_iv_curve(model, linestyle=(0, (1, 10)), label="Closest Guess") - if infeas_max > 0.5 or infeas_sum > 0.5: - all_cec_modules_df.loc[i, 'Error'] = "Infeasibility > 0.5" - if plotting: - plt.close() - continue + res, scaled_model = solve_model(model, solver, tee=False) + + if plotting: + if res and res.solver.status == 'ok': + plot_iv_curve(model, linestyle='-', label="Optimal", plot_anchors=True) + else: + plot_iv_curve(model, linestyle='dotted', label="Approx", plot_anchors=True) - a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) - all_cec_modules_df.loc[i, iv_diff_cols] = [d_I_sc, d_I_mp, d_V_mp, d_P_mp] - all_cec_modules_df.loc[i, model_param_cols] = [a, Il, Io, Rs, Rsh, Adj] - all_cec_modules_df.loc[i, 'Error'] = None + d_I_sc, d_I_mp, d_V_mp, d_P_mp = get_curve_diffs(r, model) + + infeas_sum = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).sum() + infeas_max = np.abs([d_I_sc, d_I_mp, d_V_mp, d_P_mp]).max() + if infeas_max > 0.5 or infeas_sum > 0.5: + r['Error'] = "Infeasibility > 0.5" if plotting: - plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') - plt.xlabel("Voltage") - plt.ylabel("Current") - plt.title(f"{r['Manufacturer']} {r['Model Number']}") - plt.tight_layout() - plt.savefig(plot_output_path / f"IV_curve_{i}.png") plt.close() - return all_cec_modules_df + return r + + a, Il, Io, Rs, Rsh, Adj = get_params_from_model(scaled_model) + for col, val in zip(iv_diff_cols, [d_I_sc, d_I_mp, d_V_mp, d_P_mp]): + r[col] = val + for col, val in zip(model_param_cols, [a, Il, Io, Rs, Rsh, Adj]): + r[col] = val + r['Error'] = None + + if plotting: + plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left') + plt.xlabel("Voltage") + plt.ylabel("Current") + plt.title(f"{r['Manufacturer']} {r['Model Number']}") + plt.tight_layout() + plt.savefig(plot_output_path / f"IV_curve_{i}.png") + plt.close() + return r + + +def parallel_run_solve_bootstrapping_reduced(solved_df, unsolved_df, plotting=False): + with mp.Pool(NUM_OF_WORKERS) as pool: + results = [pool.apply_async(run_solve_bootstrapping_reduced, [idx, row, solved_df, plotting]) for idx, row in unsolved_df.iterrows()] + results = [row.get() for row in results] + df = pd.DataFrame(results) + return df + +def sequential_run_solve_bootstrapping_reduced(solved_df, unsolved_df, plotting=False): + solver = pyo.SolverFactory('ipopt') + + results = [] + for i, r in unsolved_df.iterrows(): + row = run_solve_bootstrapping_reduced(i, r, solved_df, plotting, solver) + results.append(row) + return pd.DataFrame(results) def solve_bootstrapping_reduced_approx(all_cec_modules_df, solved_df, unsolved_df): @@ -848,8 +908,8 @@ def create_model_with_solution(sol_row): print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Solving for {len(all_cec_modules_df)} Modules") # solve attempt #1 - print("Starting First Pass Solve") - all_cec_modules_df = run_solve_first_pass(all_cec_modules_df) + print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting First Pass Solve") + all_cec_modules_df = parallel_run_solve_first_pass(all_cec_modules_df, plotting) all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] @@ -857,8 +917,9 @@ def create_model_with_solution(sol_row): print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Solved {len(solved_df)}. {len(unsolved_df)} remaining.") # solve with bootstrapping - print("{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Second Pass Solve") - all_cec_modules_df = run_solve_bootstrapping(all_cec_modules_df, solved_df, unsolved_df) + print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Second Pass Solve") + df = parallel_run_solve_bootstrapping(solved_df, unsolved_df, plotting) + all_cec_modules_df = pd.concat([solved_df, df]) all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] @@ -866,8 +927,9 @@ def create_model_with_solution(sol_row): print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Solved {len(solved_df)}. {len(unsolved_df)} remaining.") # solve with bootstrapping & reduced number of temperature samples - print("{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Final Pass Solve") - all_cec_modules_df = run_solve_bootstrapping_reduced(all_cec_modules_df, solved_df, unsolved_df) + print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Final Pass Solve") + df = parallel_run_solve_bootstrapping_reduced(solved_df, unsolved_df, plotting) + all_cec_modules_df = pd.concat([solved_df, df]) all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] From 2c054b6a04bafdd683cec6dd1b937d214d2ee556 Mon Sep 17 00:00:00 2001 From: Darice Date: Wed, 28 May 2025 09:10:12 -0600 Subject: [PATCH 15/17] updates from cwhanse comments --- files/SixParSolve.py | 12 +++++++++--- tests/test_SixParSolve.py | 8 ++++---- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 36028d9..b5601d1 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -25,10 +25,14 @@ model_param_cols_ssc = ["a_ssc", "Il_ssc", "Io_ssc","Rs_ssc","Rsh_ssc","Adj_ssc"] iv_diff_cols = ['d_Isc', 'd_Imp', 'd_Vmp', 'd_Pmp'] -# 6par_solve.h L 423 + def current_at_voltage_cec(Vmodule, IL_ref, IO_ref, RS, A_ref, RSH_ref, I_mp_ref): """ Solve for the current at the voltage Vmodule, given the single diode parameters which could have been modified for non-STC condition + + Used for plotting only, so numerical precision is only 1e-4. Adapted from ssc/shared/6par_solve.h L 423 + + For better stability and precision, use pvlib.pvsystem.i_from_v(v, IL_oper, IO_oper, Rs, Rsh_oper, A_oper) """ F = 0 Fprime = 0 @@ -73,8 +77,9 @@ def cec_model_params_at_condition(model, Irr, T_cell_K): A_oper = a * T_cell_K / Tc_ref EG = eg0 * (1-0.0002677*(T_cell_K-Tc_ref)) + k = pyo.value(m.k) # Boltzmann constant eV/K # instead of 1/KB, is 11600 in L129 - IO_oper = Io * np.power(T_cell_K/Tc_ref, 3) * np.exp( 11600 *(eg0/Tc_ref - EG/T_cell_K) ) + IO_oper = Io * np.power(T_cell_K/Tc_ref, 3) * np.exp(eg0 / (k*(Tc_ref)) - (EG / (k*(T_cell_K)))) Rsh_oper = Rsh*(I_ref/Irr) IL_oper = Irr/I_ref *( Il + muIsc*(T_cell_K-Tc_ref) ) @@ -170,6 +175,7 @@ def create_model(gamma_curve_dt=3): m.bVoc = pyo.Param(domain=pyo.Reals, mutable=True, units=pyo.units.V/pyo.units.K) m.gPmp = pyo.Param(domain=pyo.Reals, mutable=True, units=pyo.units.percent/pyo.units.K) m.Egref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=1.121) + m.k = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, initialize=8.617332478e-05, units=pyo.units.eV/pyo.units.K) m.Tref = pyo.Param(domain=pyo.NonNegativeReals, mutable=True, units=pyo.units.K) # Module6ParNonlinear @@ -209,7 +215,7 @@ def gamma_blocks(b, i): # PTnonlinear b.a_Tc = pyo.Expression(expr=m.par.a * b.Tc / m.Tref) b.Eg_Tc = pyo.Expression(expr=m.Egref * (1-0.0002677*(b.Tc-m.Tref))) - b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp( 11600 * (m.Egref/m.Tref - b.Eg_Tc/b.Tc))) + b.Io_Tc = pyo.Expression(expr=m.par.Io* ( b.Tc/m.Tref)**3 * pyo.exp(b.Eg_Tc / (m.k*(m.Tref)) - (m.Egref / (m.k*(b.Tc))))) b.Il_Tc = pyo.Expression(expr=m.par.Il + (m.aIsc*(1-m.par.Adj/100))*(b.Tc-m.Tref)) b.f_5 = pyo.Constraint(expr=b.Vmp_Tc *( b.Io_Tc/b.a_Tc*pyo.exp( (b.Vmp_Tc+b.Imp_Tc*m.par.Rs)/b.a_Tc ) + 1/m.par.Rsh ) / ( 1 + m.par.Rs/m.par.Rsh + b.Io_Tc*m.par.Rs/b.a_Tc*pyo.exp( (b.Vmp_Tc+b.Imp_Tc*m.par.Rs)/b.a_Tc ) ) - b.Imp_Tc == 0) diff --git a/tests/test_SixParSolve.py b/tests/test_SixParSolve.py index 07083de..9124740 100644 --- a/tests/test_SixParSolve.py +++ b/tests/test_SixParSolve.py @@ -45,11 +45,11 @@ def define_default_model(set_initial_values=True): def test_default_sam_cec_user(): model = define_default_model() - IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(model, 1000, 25+275.15) - assert IL_oper == pytest.approx(6.059759, rel=1e-5) - assert IO_oper == pytest.approx(1.1674203993060455e-10, rel=1e-5) + IL_oper, IO_oper, Rs, A_oper, Rsh_oper = cec_model_params_at_condition(model, 1000, 30+273.15) + assert IL_oper == pytest.approx(6.0683977849785, rel=1e-5) + assert IO_oper == pytest.approx(1.9115026653318136e-10, rel=1e-5) assert Rs == pytest.approx(0.3081202, rel=1e-5) - assert A_oper == pytest.approx(2.5948906, rel=1e-5) + assert A_oper == pytest.approx(2.6208265638101627, rel=1e-5) assert Rsh_oper == pytest.approx(500.069, rel=1e-5) plot_iv_curve(model) From 5fff193fb2ead952a8132a6e437b47888b0ecbcc Mon Sep 17 00:00:00 2001 From: Darice Date: Wed, 28 May 2025 09:20:10 -0600 Subject: [PATCH 16/17] sort index for output csv --- files/SixParSolve.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index b5601d1..7a5efa8 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -925,7 +925,7 @@ def create_model_with_solution(sol_row): # solve with bootstrapping print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Second Pass Solve") df = parallel_run_solve_bootstrapping(solved_df, unsolved_df, plotting) - all_cec_modules_df = pd.concat([solved_df, df]) + all_cec_modules_df = pd.concat([solved_df, df]).sort_index() all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] @@ -935,7 +935,7 @@ def create_model_with_solution(sol_row): # solve with bootstrapping & reduced number of temperature samples print(f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}: Starting Final Pass Solve") df = parallel_run_solve_bootstrapping_reduced(solved_df, unsolved_df, plotting) - all_cec_modules_df = pd.concat([solved_df, df]) + all_cec_modules_df = pd.concat([solved_df, df]).sort_index() all_cec_modules_df.to_csv(f"cec_modules_params_{filename_date}.csv", index=False) solved_df = all_cec_modules_df[~all_cec_modules_df['Rsh_py'].isna()] From 03c1a780b4db21180899f6cf017ee02a53cefb62 Mon Sep 17 00:00:00 2001 From: Darice Date: Wed, 28 May 2025 09:45:37 -0600 Subject: [PATCH 17/17] add info for IPOPT install --- files/SixParSolve.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/files/SixParSolve.py b/files/SixParSolve.py index 7a5efa8..b88bb69 100644 --- a/files/SixParSolve.py +++ b/files/SixParSolve.py @@ -892,6 +892,9 @@ def create_model_with_solution(sol_row): if __name__ == "__main__": + if not pyo.SolverFactory('ipopt').available(): + raise RuntimeError("IPOPT solver not available. Install it to your Python environment from conda: `conda install -c conda-forge ipopt`") + filename = None if len(sys.argv) > 1: