Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

clipping off-mpp operating point calculation #1178

Open
mikofski opened this issue Feb 26, 2021 · 48 comments
Open

clipping off-mpp operating point calculation #1178

mikofski opened this issue Feb 26, 2021 · 48 comments

Comments

@mikofski
Copy link
Member

Feature request
An algorithm that determines the off-mpp operating point for inverters that are clipping

Describe the solution you'd like
Assuming that if the converted AC out is greater than rated AC (at the ambient temperature), then the inverter will increase the voltage until the AC output, what is the DC output, DC voltage, and inverter efficiency loss

Describe alternatives you've considered
I scripted something really, REALLY basic, I don't suggest something as simplistic as this:

p_ac = pvlib.inverter.sandia(v_mpp, p_mpp, inverter)  # use MPP first time
voltage_test_range = (v_oc - v_mpp)/100.0  # divide excess voltage into 100 steps
v_dc = v_mpp
while p_ac > ac_rated:
    v_dc += voltage_test_range  # add a little
    i_dc = np.interp(v_dc, v, i)  # v and i are the full IV curve
    # NOTE: v must be monotonically increasing for np.interp to work!
    p_dc = v_dc * i_dc
    p_ac = pvlib.inverter.sandia(v_dc, p_dc, inverter)

But it would be better to redo this with a scalar root finder from scipy.optimize like brentq.

Also, it the AC rated power should be adjusted depending on the ambient temperature, but this could be separate.

Additional context
Sorry, I couldn't find any related issues/prs?
ICYMI: @abhisheksparikh

@cwhanse
Copy link
Member

cwhanse commented Feb 26, 2021

I agree, would be a nice feature. Somehow it should accept any of the inverter models, or, be some kind of extension within the existing functions so that a user doesn't have to choose to call an inverter function or this new function.

We could also open issues for:

  • inverter temperature derating (would be easy if there is a source for the %AC power / degree C relationship)
  • non-unity power factor operations (I think this is difficult, as it would require knowledge of how efficiency changes with power factor, which is not measured AFAIK).

@mikofski
Copy link
Member Author

may be related to #1199 ?

@kurt-rhee
Copy link
Contributor

I'd be happy to work on this in the medium term future. Is there a reference that I should work off of? Is a reference necessary for what would essentially be a point on a curve finder?

@cwhanse
Copy link
Member

cwhanse commented Dec 18, 2024

@kurt-rhee I think the request is to find the point on the IV curve where the AC output of the inverter matches a desired value. pvsystem.i_from_v (or its v_from_i counterpart) seems relevant, for the DC input. Some thought is needed to define how this capability faces the user - additional inputs for a standard inverter function? Or some new function that lives in pvsystem? I don't have an opinion here.

Also, if the inverter has multiple Arrays as inputs, and the Arrays have different conditions, there are likely to be multiple solutions. For now, I think we should assume a single Array and worry about the multiple MPPT case later.

@mikofski
Copy link
Member Author

mikofski commented Dec 18, 2024

Hi Kurt, I’m happy to jump on a call to discuss. I’m not sure if there are any references, but if not we could put this in the example gallery. I think the description already explains it well, but to clarify, the system must operate off-MPP when clipping, but finding this operating point is not trivial. It’s subject to two constraints:

  1. Ideally AC out from the inverter must equal the rated power (possibly temperature derated)
  2. The dc voltage and current must be on the system IV curve

So we probably need to iterate or use a bounded root finder like brentq. Perhaps we can just call this function?
ideal_offmpp_clipped_operation(system, ac_rating): -> (v_dc, i_dc)

@kandersolar
Copy link
Member

Another consideration: maximum/minimum voltage/current limitations of the inverter. If there is no suitable point to the right of the MPP, should the algorithm "fall back" to looking for a suitable point on the left (low voltage, high current) side?

Whether the algorithm interpolates a pre-computed I-V curve, or is given SDM parameters and calls i_from_v internally, must also be decided.

The list of considerations here is already long enough to make me wonder if we are "inventing" too much. Has anyone looked for a reference we could follow? I haven't.

@kurt-rhee
Copy link
Contributor

I also haven't been able to find a reference. Both PlantPredict and SolarFarmer have some notes about keeping the inverter input voltage and current within some operating window, but neither have references.

https://mysoftware.dnv.com/download/public/renewables/solarfarmer/manuals/latest/CalcRef/Inverter/Inverter.html

https://terabase.atlassian.net/servicedesk/customer/portal/3/article/1292009483

@mikofski
Copy link
Member Author

Good points about minimum voltage and other constraints. PVsyst also has help page with nice diagram.

I propose one of two options:

  1. Without a reference this is called the “ideal” operation. IMHO opinion this is sufficient and no reference is required because we are simply solving a set of constraints.
  2. put this in the example gallery

IMHO the reference situation is getting out of hand. Solving a set of equations is math. Only implementations of algorithms should require a reference. We are not describing a sky model here, the two situations are not comparable. There are generally available references for math and other basic solar facts.

@adriesse
Copy link
Member

Good references serve different useful purposes, but perhaps there's no need to get hung up on that right now if the fall-back is example gallery? Once there is a solution, we'll see how straight-forward it is.

@markcampanelli
Copy link
Contributor

This is interesting: Essentially we are putting a (temperature-dependent) upper bound on the quantity that we are trying to maximize. I don’t think anything a-priori prevents one from defining this like any other nonlinear constraint https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.NonlinearConstraint.html#scipy.optimize.NonlinearConstraint, and then solving for the maximum power as before, but under this additional constraint. I might try to add this feature to the MPPT-under-mismatch algorithm that I put together in #1923. Of course, the devil might be in the details!

@kurt-rhee
Copy link
Contributor

Okay finally getting around to this.

My initial thought is to create a function to solve for maximum power at the given ambient temperature #1199 and then I like @mikofski 's brentq solution for calculating the operating point given the constraint.

In a utility scale photovoltaic site with a central inverter, a given inverter will be connected to a set of combiners. Each of these combiners will then be connected to a set of modules which are not necessarily unique. For example, even intra-combiner there can be strings that have different bin classes.

This means that we will need to combine all of the IV curves together into a single IV curve before solving for the point that satisfies the maximum power constraint. This would then necessitate that a user provide a given IV curve in a 2D array of voltages and their corresponding current values.

Perhaps I am overthinking this though, looking forward to hearing everyone's thoughts.

@mikofski
Copy link
Member Author

mikofski commented Jan 7, 2025

Regarding…

This means that we will need to combine all of the IV curves together into a single IV curve before solving for the point that satisfies the maximum power constraint.

While this is true, if we only consider the “ideal” case, the situation when not all irradiance, module, & string conditions are known, could we first test the general case where the entire system is uniform? Then afterward we could provide particular solutions for more detailed scenarios? Would you be open to this approach first?

@cwhanse
Copy link
Member

cwhanse commented Jan 7, 2025

combine all of the IV curves together into a single IV curve before solving for the point that satisfies the maximum power constraint.

That's how a single MPPT works, but I don't think that's how clipping works for the whole inverter. I'm no expert by any means, but my understanding is that each MPPT applies the DC current and power limits. In contrast, the temperature limit is applied to the entire device, and I don't know how a temperature curtailment is apportioned among the MPPTs.

@mikofski
Copy link
Member Author

mikofski commented Jan 7, 2025

Maybe for simplicity sake we can descope away from temperate derating for now and just find the voltage that allows the inverter to operate within its constraints. Take this step by step and solve the simplest most idealistic case first. Then add complexity one at a time. Agile like, okay?

@wholmgren
Copy link
Member

Whether the algorithm interpolates a pre-computed I-V curve, or is given SDM parameters and calls i_from_v internally, must also be decided.

I think the function should accept an IV curve (calculated with whatever algorithms and complexity the user desires), apply the inverter parameters and logic, and output the operating point.

@adriesse
Copy link
Member

adriesse commented Jan 7, 2025

combine all of the IV curves together into a single IV curve before solving for the point that satisfies the maximum power constraint.

That's how a single MPPT works, but I don't think that's how clipping works for the whole inverter. I'm no expert by any means, but my understanding is that each MPPT applies the DC current and power limits. In contrast, the temperature limit is applied to the entire device, and I don't know how a temperature curtailment is apportioned among the MPPTs.

The total inverter power limit is probably the most frequent on to be reached and therefore the most important. Temperature de-rating is no more complicated: it just reduces the inverter total power limit value.

Typically the sum of individual MPPT current and power limits is greater than the total inverter current and power limits, so if you have multiple inputs, there has to be some logic to manage that.

@cwhanse
Copy link
Member

cwhanse commented Jan 7, 2025

Maybe for simplicity sake we can descope away from temperate derating for now and just find the voltage that allows the inverter to operate within its constraints. Take this step by step and solve the simplest most idealistic case first. Then as complexity one at a time. Agile like, okay?

Sure, but in this case it would be wise to keep DC current limits and temperature derating in mind as we alter the MPP code, since we intend to add those other capabilities.

@kurt-rhee
Copy link
Contributor

On a related note, how do people generally decide how many IV points to calculate when combining IV curves?

@mikofski
Copy link
Member Author

mikofski commented Jan 8, 2025

About 10,000 points seems reasonable. I tend to concentrate them around knees like the max power point.

@cwhanse
Copy link
Member

cwhanse commented Jan 8, 2025

10,000

That's a point every 10mV. I think one could get a satisfactory curve with 1000 points, or even fewer.

@mikofski
Copy link
Member Author

mikofski commented Jan 9, 2025

A string could be 1500[V] so 10,000 points evenly spaced is 150[mV], but between Vmp & Voc, current changes very fast. For example, a Longi LR5-72HBD-550M has only 8[V] difference between Vmp & Voc at STC, but changes by 13[A]. If an inverter has strings of 30-qty modules (16.5[kW]), then a 3.3[MW] inverter would have 200 strings, so in a space of 240[V] current changes by 2,600[A]! This region, from Vmp to Voc, will be where off-MPP operation happens, so having a lot of points helps keep current resolution more precise in my opinion.

@cwhanse
Copy link
Member

cwhanse commented Jan 9, 2025

Good points, I wasn't thinking about string-level IV curves.

@kurt-rhee
Copy link
Contributor

This is great insight. @mikofski is there a good algorithm for concentrating iv points around the MPPs?

We can reduce the amount of points needed by only considering voltages within the inverter operating window (excluding 0V to Vmin), but the above comment from mark definitely helps me understand what sort of bin? size we need between iv points

@mikofski
Copy link
Member Author

mikofski commented Jan 9, 2025

is there a good algorithm for concentrating iv points

Unfortunately I don’t have good advice. In the past I’ve used two log-spaces reflected at Vmp but it’s not perfect. (See bishop88.) Also you may encounter many local minima, so an adaptive mesh that is directly proportional to dI/dV might work better. It’s sort of an implicit problem - you need the mesh to know what mesh to use. :/

@adriesse
Copy link
Member

adriesse commented Jan 9, 2025

In my experience the points accumulate in the right places as I build up the composite IV curve, so I never thought much about how many points in total.

@kurt-rhee
Copy link
Contributor

Forgive my ignorance, but I am not sure how to create a composite IV curve without selecting a bunch of reference voltages to calculate current at.

A part of me thinks that it might be possible to only combine points on curves that correspond with another curve's MPP, but I am not sure.

Appreciate any insights!

@cwhanse
Copy link
Member

cwhanse commented Jan 10, 2025

To "add" IV curves for devices in series, I think in terms of drawing horizontal lines, and adding the voltages at each current value. That requires extending each IV curve to currents up to the maximum Isc for all curves. A model for this negative bias part of the curve is thus needed.

The bishop88 functions provide a model. But there are convergence problems when solving the extended IV curve using the Newton method (#1757), and the BrentQ method isn't vectorized, so calculations would be very slow for a large set of curves.

I've taken to simply linearly extending the IV curve to the breakdown voltage and then drawing a vertical asymptote. That simplified calculating current at negative voltages and makes the addition of many curves (or time series of many curves) feasible, without a meaningful loss in precision.

#1781 is a start at a function to add IV curves. It got hung up on the convergence and performance issues with the bishop88 functions. It would be great to see that pushed forward.

@wholmgren
Copy link
Member

It seems to me that the "clipping off-mpp operation point calculation" function should be tested against a handful of IV curves ranging from simple to complex and perhaps with a range of discretization. But if the function merely accepts and IV curve as an input, which I believe is the right way forward, then for now these IV curve calculation problems can be solved by hacky one off code instead of robust pvlib implementations.

@mikofski
Copy link
Member Author

Is adding IV curves in scope for this PR? I was thinking the function would only take in a single IV curve then find the operating point that meets the inverters constraints (if not @ MPP) - maybe this is oversimplifying but even multi MPPT inverters have to first handle this individually for each channel right?

@adriesse
Copy link
Member

On a related note, how do people generally decide how many IV points to calculate when combining IV curves?

@kurt-rhee perhaps you can start a separate discussion on this topic, if needed.

@markcampanelli
Copy link
Contributor

markcampanelli commented Jan 11, 2025

Here is a DC-side solution for a single 60-cell PV module that uses scipy.optimize.minimize's method=''trust-constr'. Lot's of rough edges here (edge cases, feasibility assurance for the initial condition, tolerance tuning, performance, etc.), but it should get the idea across. It runs pretty quick on my old laptop. I would hope that this approach could be usefully extended to the more complicated case where there is coupling between the AC inverter efficiency and the DC input, which I presume complicates the simple max_power_derate used here.

One might wonder why I choose to use current as the variable to optimize over, which leads to a nonlinear constraint on voltage. The reason is because I'd like to extend this to the case where panels in a string are mismatched.

import numpy
import pvlib.pvsystem
import scipy.constants
import scipy.optimize

photocurrent = 6.2
saturation_current = 1.0e-8
resistance_series = 0.0001 
resistance_shunt = 5000
n = 1.1
n_s = 60
k_B_J_per_K = scipy.constants.value("Boltzmann constant")
T_K = scipy.constants.convert_temperature(25.0, "Celsius", "Kelvin")
q_C = scipy.constants.value("elementary charge")
nNsVth = n * n_s * k_B_J_per_K * T_K / q_C

args = (
    photocurrent,
    saturation_current,
    resistance_series,
    resistance_shunt,
    nNsVth,
)

max_power_point_unconstrained = pvlib.pvsystem.max_power_point(*args)
max_power_derate = 0.95
p_mp_max = max_power_derate*max_power_point_unconstrained["p_mp"]
v_mp_min = max_power_point_unconstrained["v_mp"]
i_mp_ic = 0.9*max_power_derate*max_power_point_unconstrained["i_mp"]  # Must be feasible.


# Note closure over args (could try a functools.partial).
def negative_power_partial(current):
    """
    Compute negative of power generated by device at specified current.
    Designed for use by scipy.optimize.minimize.
    """
    return -current * pvlib.pvsystem.v_from_i(current, *args)


# Note closure over args (could try a functools.partial).
def voltage_at_current_partial(current):
    return pvlib.pvsystem.v_from_i(current, *args)


nonlinear_constraint_power = scipy.optimize.NonlinearConstraint(
    negative_power_partial,
    -p_mp_max,
    0,
    jac='3-point',
    keep_feasible=True,  # Need ensure ic is feasible if we make this True.
)

nonlinear_constraint_voltage = scipy.optimize.NonlinearConstraint(
    voltage_at_current_partial,
    v_mp_min,
    numpy.inf,
    jac='3-point',
    keep_feasible=True,  # Need ensure ic is feasible if we make this True.
)

sol = scipy.optimize.minimize(
    negative_power_partial,
    i_mp_ic,
    jac='3-point',
    method='trust-constr',
    constraints=(nonlinear_constraint_power, nonlinear_constraint_voltage),
)

print(
    "Maximum power point unconstrained:\n"
    f"i_mp={max_power_point_unconstrained['i_mp']}, "
    f"p_mp={max_power_point_unconstrained['p_mp']}, "
    f"v_mp={max_power_point_unconstrained['v_mp']}."
)

print("")

print(
    f"Maximum power point constrained to p_mp <= {p_mp_max} and v_mp >= {v_mp_min} "
    f"(max_power_derate={max_power_derate}):\n"
    f"i_mp={sol.x[0]}, "
    f"p_mp={-sol.fun}, "
    f"v_mp={-sol.fun/sol.x[0]}."
)

Should output something like:

Maximum power point unconstrained:
i_mp=5.856590872283306, p_mp=172.14769430100824, v_mp=29.39383987290427.

Maximum power point constrained to p_mp <= 163.54030958595783 and v_mp >= 29.39383987290427 (max_power_derate=0.95):
i_mp=5.252975152617371, p_mp=163.53629540662442, v_mp=31.13212810937057.

@markcampanelli
Copy link
Contributor

markcampanelli commented Jan 11, 2025

Here's a version that directly incorporates an inverter's AC power limit. This code takes a bit longer to run (~8s) and may need to be optimized to be practical.

What's pretty fun is that I think one can incorporate more practical nonlinear constraints at will, such as the MPP DC voltage window, and maximum DC and/or AC current limits.

import numpy
import pvlib.pvsystem
import scipy.constants
import scipy.optimize

# module
photocurrent = 7.0
saturation_current = 1.0e-8
resistance_series = 0.0001 
resistance_shunt = 5000
n = 1.1
n_s = 60
k_B_J_per_K = scipy.constants.value("Boltzmann constant")
T_K = scipy.constants.convert_temperature(25.0, "Celsius", "Kelvin")
q_C = scipy.constants.value("elementary charge")
nNsVth = n * n_s * k_B_J_per_K * T_K / q_C

module_params = (
    photocurrent,
    saturation_current,
    resistance_series,
    resistance_shunt,
    nNsVth,
)

# inverter
db = pvlib.pvsystem.retrieve_sam('cecinverter')
inverter = db[db.loc["Pdco"].idxmin()]  # Get smallest-capacity inverter.

max_power_point_unconstrained = pvlib.pvsystem.max_power_point(*module_params)
v_mp_min = max_power_point_unconstrained["v_mp"]
i_mp_ic = 0.9*max_power_point_unconstrained["i_mp"]  # Must be feasible.


# Note closure over args (could try a functools.partial).
def opposite_dc_power(current):
    """
    Compute opposite of power generated by device at specified current.
    Designed for use by scipy.optimize.minimize.
    """
    return -current * pvlib.pvsystem.v_from_i(current, *module_params)


# Note closure over args (could try a functools.partial).
def ac_power_at_dc_current(current):
    """Designed for use by scipy.optimize.NonlinearConstraint."""
    v_dc = pvlib.pvsystem.v_from_i(current, *module_params)
    return pvlib.inverter.sandia(v_dc, v_dc*current, inverter)


# Note closure over args (could try a functools.partial).
def dc_voltage_at_dc_current(current):
    """Designed for use by scipy.optimize.NonlinearConstraint."""
    return pvlib.pvsystem.v_from_i(current, *module_params)


nonlinear_constraint_ac_power = scipy.optimize.NonlinearConstraint(
    ac_power_at_dc_current,
    -numpy.inf,
    inverter["Paco"],  # Could temperature-derate this.
    jac='3-point',
    keep_feasible=True,  # Need to ensure IC is feasible if we make this True.
)

nonlinear_constraint_dc_voltage = scipy.optimize.NonlinearConstraint(
    dc_voltage_at_dc_current,
    v_mp_min,
    numpy.inf,
    jac='3-point',
    keep_feasible=True,  # Need to ensure IC is feasible if we make this True.
)

sol = scipy.optimize.minimize(
    opposite_dc_power,
    i_mp_ic,
    jac='3-point',
    method='trust-constr',
    constraints=(nonlinear_constraint_ac_power, nonlinear_constraint_dc_voltage),
)

print(
    "Maximum power point unconstrained:\n"
    f"i_mp={max_power_point_unconstrained['i_mp']}, "
    f"p_mp={max_power_point_unconstrained['p_mp']}, "
    f"v_mp={max_power_point_unconstrained['v_mp']}."
)

print("")

print(
    f"Maximum power point constrained to p_ac <= {inverter['Paco']} and v_mp >= {v_mp_min}:\n"
    f"i_mp={sol.x[0]}, "
    f"p_mp={-sol.fun}, "
    f"v_mp={-sol.fun/sol.x[0]}, "
    f"p_ac={ac_power_at_dc_current(sol.x[0])}."
)

Should output something like:

[/home/campy/Documents/Employment/IMS/solar/PVfit/venv/pvlib/lib/python3.11/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py:81](https://file+.vscode-resource.vscode-cdn.net/home/campy/Documents/Employment/IMS/solar/PVfit/venv/pvlib/lib/python3.11/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py:81): UserWarning: Singular Jacobian matrix. Using SVD decomposition to perform the factorizations.
  Z, LS, Y = projections(A, factorization_method)
[/home/campy/Documents/Employment/IMS/solar/PVfit/venv/pvlib/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:551](https://file+.vscode-resource.vscode-cdn.net/home/campy/Documents/Employment/IMS/solar/PVfit/venv/pvlib/lib/python3.11/site-packages/scipy/optimize/_differentiable_functions.py:551): UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.
  self.H.update(delta_x, delta_g)
[/home/campy/Documents/Employment/IMS/solar/PVfit/venv/pvlib/lib/python3.11/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py:166](https://file+.vscode-resource.vscode-cdn.net/home/campy/Documents/Employment/IMS/solar/PVfit/venv/pvlib/lib/python3.11/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py:166): RuntimeWarning: invalid value encountered in scalar subtract
  actual_reduction = merit_function - merit_function_next
Maximum power point unconstrained:
i_mp=6.615301306318015, p_mp=195.74202445950908, v_mp=29.589283298791475.

Maximum power point constrained to p_ac <= 175.0 and v_mp >= 29.589283298791475:
i_mp=5.953771175686214, p_mp=186.36492485609355, v_mp=31.301996559283904, p_ac=175.0.

@cwhanse
Copy link
Member

cwhanse commented Jan 13, 2025

@markcampanelli do you know if that procedure be vectorized over IV curves? From what I see in the scipy documentation, I think not. The function to be minimized has to return float, I think.

@adriesse
Copy link
Member

Feature request An algorithm that determines the off-mpp operating point for inverters that are clipping

Describe the solution you'd like Assuming that if the converted AC out is greater than rated AC (at the ambient temperature), then the inverter will increase the voltage until the AC output, what is the DC output, DC voltage, and inverter efficiency loss

One thing that's missing from this feature request is the motivation. How important is it to know this voltage?

@cwhanse
Copy link
Member

cwhanse commented Jan 13, 2025

@markcampanelli's code sample gave me some inspiration in how to think about the various limits that inverters apply. There are limits on DC current and voltage, and on AC current which are usually thought of as limits on temperature and AC power.

Limits on DC current and DC voltage:

  • DC current has a maximum. When array current exceeds this maximum, inverters shift to higher voltage.
  • inverters have a start-up or turn-on voltage. Below this limit, the inverter presents an open circuit to the array.
  • MPPT voltage has a range. Inverters can (and do) operate outside this range, but the conversion efficiency is not routinely (ever?) characterized, and we have no applicable model outside the MPPT range.

Importantly, DC constraints apply independently to each MPPT.

Limits on AC current:

  • The usual AC power limit is (I believe) actually a limit on AC current.
  • The temperature limit also appears to be applied by constraining AC current, see e.g. SolarEdge device.

The AC power limit applies to the whole inverter. It is not clear how a power derate is apportioned among several MPPTs.

All of these limits are likely implemented in the inverter controls, rather than in hardware, because the response to each is to shift the DC voltage. It is not clear if the constraints are applied in sequence or simultaneously. But I don't think that matters for performance modeling; we're interested in the voltage-current-power after all these constraints are applied.

The DC current and voltage limits are simple constraints on the domain for voltage and current, which define a segment of the IV curve. The AC power limit affects which point is selected within the IV curve segment.

Because there is no inverter temperature model (AFAIK), and the relationship between inverter temperature and output AC current is unknown, this constraint should be set aside for now. But if the pvlib function that determines the operating point can accept a power limit constraint, adding a constraint on AC current for temperature shouldn't be difficult.

@markcampanelli
Copy link
Contributor

markcampanelli commented Jan 13, 2025

@cwhanse I did not look into vectorization, but I would not be surprised if you are correct that this more complicated solver is not vectorized "out of the box". I suppose this was part of the "performance" items I said were left to investigate.

I agree that there are opportunities here to implement several additional inverter constraints, and I'm excited to see what you can come up with :)! My thoughts were along the lines of some sort of a higher-level "mppt string-input computation" function that:

  1. Takes as inputs a module's configured v_from_i function (equiv. i_from_v?) and an inverter's configured dc_to_ac function, and the number of cells in the string,
  2. Reads the config for the inverter and constructs as many linear/nonlinear constraints as appropriate for the provided inverter info,
  3. Solves the constrained DC power maximization, outputting both the DC voltage and current for the string and AC power (and ???).

In this fashion, we would be setting up an adaptable interface that (hopefully) would accommodate a "mix and match" of module IV and inverter DC-to-AC models to compute string inputs to MPPTs. I think this might even be fairly easy to extend to a mismatched case, as I sketched out in #1923.

Also, one tricky case (which may be a rounding error in the big picture) is the "turn on DC voltage" for the MPPT that may be different than the lower end of the voltage window. I cannot think of how to catch this other than with a post-processing filter that looks at the entire time-series of voltages. More generally, there might need to be a filtering step that lives AFTER the power-maximization function.

Lastly, sorry for another big code listing here, but I simplified my above code a bit to show how some nonlinear constraints can be combined (see esp. ac_power_and_dc_voltage_at_dc_current). However, this approach might be harder to "mix and match" than just providing a tuple of individual constraints, and I did not do a performance comparison.

import numpy
import pvlib.pvsystem
import scipy.constants
import scipy.optimize

# module
photocurrent = 7.0
saturation_current = 1.0e-8
resistance_series = 0.0001 
resistance_shunt = 5000
n = 1.1
n_s = 60
k_B_J_per_K = scipy.constants.value("Boltzmann constant")
T_K = scipy.constants.convert_temperature(25.0, "Celsius", "Kelvin")
q_C = scipy.constants.value("elementary charge")
nNsVth = n * n_s * k_B_J_per_K * T_K / q_C

module_params = (
    photocurrent,
    saturation_current,
    resistance_series,
    resistance_shunt,
    nNsVth,
)

# inverter
db = pvlib.pvsystem.retrieve_sam('cecinverter')
inverter = db[db.loc["Pdco"].idxmin()]  # Get smallest-capacity inverter.

max_power_point_unconstrained = pvlib.pvsystem.max_power_point(*module_params)
v_mp_min = max_power_point_unconstrained["v_mp"]
i_mp_ic = 0.9*max_power_point_unconstrained["i_mp"]  # Must be feasible.


# Note closure over args (could try a functools.partial).
def opposite_dc_power(current):
    """
    Compute opposite of power generated by device at specified DC current.

    Designed for use by scipy.optimize.minimize.
    """
    return -current * pvlib.pvsystem.v_from_i(current, *module_params)


# Note closure over args (could try a functools.partial).
def ac_power_and_dc_voltage_at_dc_current(current):
    """
    Compute inverter-constrained quantities at specified DC current.

    Designed for use by scipy.optimize.NonlinearConstraint.
    """
    v_dc = pvlib.pvsystem.v_from_i(current, *module_params)

    return numpy.concat((pvlib.inverter.sandia(v_dc, v_dc*current, inverter), v_dc))


nonlinear_constraint = scipy.optimize.NonlinearConstraint(
    ac_power_and_dc_voltage_at_dc_current,
    (0, v_mp_min),
    (inverter["Paco"], numpy.inf),  # Could temperature-derate Paco.
    jac='3-point',
    keep_feasible=True,  # Need to ensure IC is feasible if we make this True.
)

sol = scipy.optimize.minimize(
    opposite_dc_power,
    i_mp_ic,
    jac='3-point',
    method='trust-constr',
    constraints=nonlinear_constraint,
)

print(
    "Maximum power point unconstrained:\n"
    f"i_mp={max_power_point_unconstrained['i_mp']}, "
    f"p_mp={max_power_point_unconstrained['p_mp']}, "
    f"v_mp={max_power_point_unconstrained['v_mp']}."
)

print("")

print(
    f"Maximum power point constrained to p_ac <= {inverter['Paco']} and v_mp >= {v_mp_min}:\n"
    f"i_mp={sol.x[0]}, "
    f"p_mp={-sol.fun}, "
    f"v_mp={-sol.fun/sol.x[0]}, "
    f"p_ac={pvlib.inverter.sandia(-sol.fun/sol.x[0], -sol.fun/sol.x[0]*sol.x[0], inverter)}."
)

@mikofski
Copy link
Member Author

This could be a chance to use SciPy Cython Optimize toolbox. Performance can be 30x native map() when nogil is used. See https://gist.github.com/mikofski/ac30065073d0d32d6ea3569f6e24e5ec

@cwhanse
Copy link
Member

cwhanse commented Jan 14, 2025

Continuation of my comment about constraints: finding the maximum power point subject to all three types of constraint (DC voltage, DC current and AC power) is, I think, a convex optimization, when the IV curve does not have mismatch features.

For one MPPT, finding the maximum DC power point is a convex optimization problem: the downward concave shape of an IV curve is convex, and the feasible set of solutions (defined by the above linear constraints) is also convex.

The inverter operating point is not synonymous with the maximum power point: the difference is, the inverter operating point is constrained by limits on DC voltage, DC current, AC current and AC power. So shift terminology here.

The standard form of the optimization problem for the inverter operating point, with DC constraints, is:

Minimize f(Vdc) = - I(Vdc) * Vdc

subject to:
Vdc <= Vmax
-Vdc <= -Vmin
I(Vdc) < ImaxDC

For one MPPT, finding the inverter operating point that also respects the AC power limit may be convex if the inverter model satisfies certain conditions (convex and monotone). I think the Sandia model satisfies these conditions for most parameter sets, but this would need confirmation.

If so, the operating point problem in standard form is

Minimize g(Vdc) = -inverter(I(Vdc), Vdc))

subject to:

Vdc <= Vmax
-Vdc <= -Vmin
I(Vdc) < ImaxDC
g(Vdc) <= PmaxAC

For multiple MPPTs, we need an assumption for how an inverter would allocate AC power reductions among the MPPTs. Because power is additive over MPPTs, I still have hope that the optimization problem will be convex.

Because each timestep is independently solved, vectorization is only a matter of adding rows to a design matrix.

Why focus on convexity? Because convex optimization is an extensively researched topic, and large-scale optimizers are available. Is there a optimizer that is free to use and easy to integrate with pvlib? I don't know, yet.

@adriesse
Copy link
Member

Continuation of my comment about constraints: finding the maximum power point subject to all three types of constraint (DC voltage, DC current and AC power) is, I think, a convex optimization, when the IV curve does not have mismatch features.

Why focus on convexity? Because convex optimization is an extensively researched topic, and large-scale optimizers are available.

Just wondering, are you thinking you may not need to do the optimization when the IV curve has mismatch features, or just putting that challenge off until later?

@cwhanse
Copy link
Member

cwhanse commented Jan 14, 2025

putting that challenge off until later?

This. It may be that the Vdc domain could be divided into subsets within which the problem is convex, or that may prove impractical.

@cwhanse
Copy link
Member

cwhanse commented Jan 15, 2025

There are tools in SciPy optimize for convex optimization, for example brentq & newton & toms748 & chandrupatla see:

Those are root finders and minimizers for scalar functions, which don't leverage the convexity of the domain or the objective function. Only the Newton method has been natively vectorized (as you know!)

@markcampanelli
Copy link
Contributor

I haven't fully digested the recent comments here, but when it comes to an extension to MPPT tracking, it seems that we could start to wade into a situation where we often do not know the "rules" of how any particular inverter will actually track a mismatched string as the weather evolves on any given day. So, the question becomes: Can we come up with a reasonable approximation that "overall" does better than the matched assumption? Indeed, I view the apparent fact that inverters generally make DC voltage larger in order to limit AC voltage one of these "rules".

@cwhanse
Copy link
Member

cwhanse commented Jan 15, 2025

What do you mean by the "matched assumption"? That series-connected devices conduct a common current?

@mikofski
Copy link
Member Author

mikofski commented Jan 16, 2025

fsolve is for convex problems

https://docs.scipy.org/doc/scipy-1.15.0/reference/generated/scipy.optimize.fsolve.html

Under the hood it calls the FORTRAN minpack library hybrdj function

@cwhanse
Copy link
Member

cwhanse commented Jan 16, 2025

fsolve is for convex problems

fsolve doesn't have bounds constraints, as far as I can tell.

I think we may mean different things by "convex." fsolve finds the root of a set of equations, i.e., F(x) = 0 where x \in R^n. In concept, folsve could be applied to G', where G = F^2 and in this instance, G is convex. Solving that requires 1) F to be differentiable and 2) x to be in the interior of the feasible set. The second condition is not likely to hold for finding the operational point for an inverter, which will sometimes be found on the boundary of the feasible set, e.g., when Vdc = Vmppt,max.

A convex optimization problem has both a convex objective (loss) function and a convex feasible set of solutions.

@adriesse
Copy link
Member

Actually, the feasible set of solutions may not be convex with a maximum power constraint because there can be a feasible solution at a high voltage and low current as well as at low voltage and high current.

@kurt-rhee
Copy link
Contributor

Status update for everyone, I got asked to work on a different project so I am going to have to put this down for now. Will hopefully return in the future.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

7 participants