Price taker model for DISPATCHES, Rehashed

[djlaky]

Fixes

Compared to #1201, operational constraints mathematical form was corrected. Unnecessary functions were removed/merged. Additional user flexibility was added for constructing cost objectives.

Summary/Motivation:

Resurrecting #1201 to finish price taker framework in accordance with project milestones.

Framework allows the user to construct price-taker models for design and/or operational optimization considering time-varying market price data.

Legal Acknowledgement

By contributing to this software project, I agree to the following terms and conditions for my contribution:

1. I agree my contributions are submitted under the license terms described in the LICENSE.txt file at the top level of this directory.
2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

[andrewlee94]

Just to check, what are these dependencies required for and are they absolutely necessary?

[djlaky]

They're used for KMeans clustering to reduce a large set of days (~365 days) to a smaller number (output an 'optimal' clusters function, or user specifiedvalue):

idaes-pse/idaes/apps/grid_integration/multiperiod/price_taker_model.py

Lines 170 to 178 in 652d1cc

	for k in k_values:
	kmeans = KMeans(n_clusters=k).fit(daily_data.transpose())
	inertia_values.append(kmeans.inertia_)
	# Identify the "elbow point"
	elbow_point = KneeLocator(
	k_values, inertia_values, curve="convex", direction="decreasing"
	)
	n_clusters = elbow_point.knee

and:

idaes-pse/idaes/apps/grid_integration/multiperiod/price_taker_model.py

Line 239 in 652d1cc

kmeans = KMeans(n_clusters=n_clusters).fit(daily_data.transpose())

@radhakrishnatg, Can you comment on if these dependencies (kneed, scikit-learn) are necessary?

[andrewlee94]

This is not something for this PR, but one thing I have wanted to comment on for a while is to make sure that the DISPATCHES team is aware the IDAES steady-state flowsheets support time indexing.

m.fs = FlowsheetBlock(dynamic=False, time_set=[0, 1, 2, ...])

I do not think this tool has been using this capability, and it might help simplify things (or conversely make things more complicated). However, this would allow you to e.g. initialize all states (StateBlocks) in the multi-period model in one go using the existing initialization frameworks.

[djlaky]

We will note this for our continued development, thank you.

[radhakrishnatg]

@andrewlee94 Thank you for the note! I was not aware of time indexing. In DISPATCHES, we create one instance of the flowsheet, initialize it and clone it. Also, in most cases, models tend be linear/quadratic, so initialization is not needed. If we want to use the time index, then I believe we have to make a lot of changes to the MultiPeriodModel class. Maybe for another PR.

[andrewlee94]

Definitely another PR.

For reference, state block initialization will handle the indexing just fine and will initialize all time points in one go for you. @Robbybp also has some dynamic initialization tools that copy results from one state to another that might work here as well.

[Robbybp]

Yes, I have some tools in pyomo-mpc for working with dynamic models. These should be useful for working with time-indexed steady state flowsheets, as this is very close to the use case they were designed for. To my basic understanding, MultiPeriodModel and the related tooling in Dispatches seems like they are solving a similar problem from a different angle.

[andrewlee94]

We generally discourage print statements inside code; using a logger is generally better as it gives control over where the output goes.

[andrewlee94]

Should this be moved into a separate method so that it can be called by the user separately?

[andrewlee94]

Looking more closely, I would recommend breaking this into a separate method. This will make testing easier as well.

[andrewlee94]

Does this really test much of the actual code? It runs the plotting code, but all this assert does is make sure the resulting plot got closed which doesn't really test the plot. Note that testing plotting code is generally very hard - we generally only test the data used to generate the plot and leave the actual plotting untested.

[djlaky]

You're correct this just tests if the plot was successfully called; it does not check if the "correct" plot or "correct" data was displayed.

I'm inclined to remove the plotting test. Please see my other review comment for further information on the optimal clusters code snippet (link).

[djlaky]

It would appear that changing platforms where the code is run or changing the version of scikit-learn can generate different results for the optimal clustering, changing the result of this plot as well.

That's why I made the test a range instead of a single value in another test:

idaes-pse/idaes/apps/grid_integration/multiperiod/tests/test_price_taker.py

Lines 76 to 88 in 652d1cc

	def test_determine_optimal_num_clusters(excel_data):
	# Added a range for optimal cluster values based on how the
	# plot appears visually. Test can be removed in the future if
	# failure occurs. This may depend on scikit-learn and kneed and
	# the interaction thereof.
	# Older versions get n_clusters = 10, Newer versions n_clusters = 11
	m = PriceTakerModel()
	daily_data = m.generate_daily_data(excel_data["BaseCaseTax"])
	n_clusters, inertia_values = m.get_optimal_n_clusters(daily_data)
	assert 9 <= n_clusters <= 15

Ultimately, whether the above test stays hinges on whether the finding of "optimal" clusters stays. I'm also inclined to keep the optimal clustering functionality.

@radhakrishnatg, can you please give us some guidance?

[andrewlee94]

This would imply that there should be bounds on the version of sklearn we support to ensure expected behaviour (or at least for testing purposes).

[radhakrishnatg]

@djlaky Not sure what this is trying to do. From what I understand, you are trying to set the LMP values which are less than 1e-4 to zero to avoid numerical issues. This is good, and we need this. But, it cannot be value <=1e-4. LMP values can also be negative, and this would set those LMPs to zero as well. We need to set LMPs to zero if their absolute value is < 1e-4.

[djlaky]

I made the change. Good catch

[radhakrishnatg]

You don't have to change it, but you can write it this way:
min_capacity * op_mode[t] <= var_commodity[t]

[andrewlee94]

I think you should make that change - for one it simplifies the code and takes it easier to read, and secondly it removes one unnecessary method call that could potentially break in the future.

[djlaky]

This was originally in place to handle cases where multi-year, or multi-day (or both) were used to specify the multiperiod model. In this case, self.mp_model.set_period is a data structure full of tuples i.e.:

instead of:
set_period = [period_1, period_2, ..., period_n]

we have:
set_period = [(period_1, year_1), (period_2, year_1), ... , (period_n, year_1), (period_1, year_2), ..., ..., (period_n, year_n)]

That's why the structure self.range_time_steps was created. In this way, you can use the .at(t) functionality to get the tuple within the self.mp_model.set_period data structure indexed by the self.range_time_steps RangeSet.

Three realistic options:

1. Leave it the way it is.
2. Make a call t = self.mp_model.set_period.at(t) at the beginning of these functions to make it more readable.
3. Update it to be more readable for now and inevitably change it back when we use representative days or multi-year simulations.

Since the multiperiod class is the way it is, I think this probably has to stay the way it is for the moment. We can address this in a future release when we fully support representative days (which will also use the tuple structure from the multiperiod class). Up to suggestions, though.

[andrewlee94]

Where and how is this rule used? My inclination would be handle this at the component generation step, and only used the indexing elements you want. That way, this rule will only ever get a single value for t.

[andrewlee94]

Also, why not have two separate indexing sets in this case, one for period and one for year? That way they would be clearly distinguished (and year can be a single element if you don't need it).

Whilst dealing with arbitrary length indexing elements is possible, it does add some complexity and edge cases you need to be able to handle. This brings me to another question, do you have tests for both cases (and any edge cases you can envision occurring)?

[radhakrishnatg]

@djlaky We have declared:

op_mode = {
            t: self.mp_model.period[t].find_component(op_blk + ".op_mode")
            for t in self.mp_model.period
        }

so, t = self.mp_model.set_period.at(t). If elements of set_period are tuples, then t is also a tuple. Right now, we only plan to support representative days. The multiyear case can be handled by creating an instance of the PriceTakerModel for each year. If single index vs. tuple is the issue, then you can construct the constraint in this manner.

m.capacity_low_lim = Constraint(m.set_period)
for t in m.set_period:
    m.capacity_low_lim.body(min_capacity * op_mode[t] <= var_commodity[t])

This can handle both single index and double index. Double check the method name once again. It has to be either body or expr or something like that.
Alternatively, you can separate it out depending on full year vs. representative days.

[radhakrishnatg]

Apparently, there is another approach:

@m.Constraint(m.set_period)
def capacity_low_lim(blk, *t):
    return min_capacity * op_mode[t] <= var_commodity[t]

This works for both single index and tuple.

[djlaky]

For constraints that have no skipping, sure, but for the ramping constraints and startup/shutdown there are issues:

https://github.com/djlaky/idaes-pse/blob/price-taker-model/idaes/apps/grid_integration/multiperiod/price_taker_model.py#L644-L656

    # The linearized ramping constraints
    @blk.Constraint(self.range_time_steps)
    def ramp_up_con(b, t):
        if t == 1:
            return Constraint.Skip
        else:
            return (
                var_ramping[self.mp_model.set_period.at(t)]
                - var_ramping[self.mp_model.set_period.at(t - 1)]
                <= (startup_rate - ramp_up_rate)
                * act_startup_rate[self.mp_model.set_period.at(t)]
                + ramp_up_rate * act_op_mode_rate[self.mp_model.set_period.at(t)]
            )

https://github.com/djlaky/idaes-pse/blob/f10ea59ad3f2a430f749748b42be8419e34982a3/idaes/apps/grid_integration/multiperiod/price_taker_model.py#L751-L759

    def Binary_relationhsip_con(b, t):
        if t == 1 or t > number_time_steps:
            return Constraint.Skip
        return (
            op_mode[self.mp_model.set_period.at(t)]
            - op_mode[self.mp_model.set_period.at(t - 1)]
            == start_up[self.mp_model.set_period.at(t)]
            - shut_down[self.mp_model.set_period.at(t)]
        )

[radhakrishnatg]

Looks good to me.

[andrewlee94]

Can you comment on why this might not exist at this point? It seems to me that this is something that should exist before you get here.

[andrewlee94]

1. We will need documentation for all the new classes - at a minimum some technical references for the API (autodocs), but and explanation with a simple demonstration of how to use the tools would be nice if possible (looking at the code, this will be very important).
2. I do not see a test file for design_and_operation_model.py. This might be rolled into another test, but I recommend having a separate test file for this to help with future maintenance.
3. In general, I think you could (and should) do a lot more unit testing of methods in isolation. From what I can see in the tests, there are not a lot of tests that target individual methods. E.g. I see no test that targets add_startup_shutdown and actually tests that it does what it is supposed to do (I see one test for a failure case, but that is it).
4. Most of the tests I see are based on what I assume is an actual case study. I would suggest thinking about whether you can construct simpler dummy cases for unit testing all the methods in isolation, and where possible mock up any necessary starting information so that you have fewer moving parts and simpler cases to test with. The full case study is good for an end-to-end component/integration test of the system, but often make it harder to do unit testing.

[andrewlee94]

What is the purpose of this test? As far as I can tell, this test would only fail if PriceTakerModel somehow prevents users from creating new attributes. If you expect seed to be created by PriceTakeModel, you should use either assert isinstace(m.seed, expected_type) or assert hasattr(m, "seed").

[andrewlee94]

Seeing as seed is a property with a custom settr, I would suggest testing for the presence of _seed and the failure cases for seed here.

[andrewlee94]

I would not use an f-string here - check for the expected string output to make sure nothing else is going wrong. Same for the above and below.

[andrewlee94]

This is one very long test involving a lot of with pytest.raises() checks. I would suggest breaking this in to multiple tests which check each possible exception individually. This will help future maintenance, as you will be able to see all exceptions in one pass; at the moment the test will fail at the first exception it encounters and other exceptions might get hidden until you fix the first. This applies to some other tests I've seen too.

[andrewlee94]

I would suggest moving the with pytest.rasies() to be only around the final piece of code you expect to fail. This way, any failures in the set up code will cause different errors to be emitted which might help with debugging.

Related to this:

1. Are there any unit tests for each of the lines of code in the set up part, to make sure these are tested in isolation? At the moment, this tests depends on a lot of moving parts and I am not sure if they have been tested sufficiently to separate a failure in the setup for the failure you are actually trying to test. This goes for a lot of the other tests here too.
2. Is there a way to dummy up the first part to avoid possible failures in the set up code?

[andrewlee94]

Pylint will probably complain about import ordering here. You need to group imports, and start with Python native funcitons, then thrid party imports, then Pyomo and finally IDAES.

[andrewlee94]

This is probably OK, but would it be better to use the IDAES logger instead?

[andrewlee94]

This would imply that there should be bounds on the version of sklearn we support to ensure expected behaviour (or at least for testing purposes).

[andrewlee94]

Looking more closely, I would recommend breaking this into a separate method. This will make testing easier as well.

[andrewlee94]

A minor comment, but I would generally put the simple checks first before trying to open the file. These checks are a lot cheaper than opening a file, and I prefer to fail sooner rather than later.

[andrewlee94]

This is not a great test (I assume this is what you were referring to with different platforms/versions giving different results). This also does not test the inertia values which would be good to include (if possible). Would it be possible to use a much simpler dummy data set to test this which gives a well defined solution for testing?

Also - this is where break the plotting function out as a separate method would help as you could have a separate test for it.

[ksbeattie]

@djlaky we're hoping to cut a release at the end of the month - next week. Do you think this will be done in time for that?

[adowling2]

We will not have time to incorporate the feedback from Andrew. If we are okay with merging without those suggestions (maybe open an issue to track proposed changes), then yes, this can probably get merged.

[ksbeattie]

Pushing this to the Aug release

[ksbeattie]

@djlaky any update on this?