Problems with unit ramp limits

[Mastomaki]

Description

Ramping down will be difficult when time resolution changes.

How to reproduce

Set the following unit and unit-output node parameters for a simple generator unit:

fix_units_on: 1
unit_capacity: 100
minimum_operating_point: 0.4
ramp_up_limit: 0.1
ramp_down_limit: 0.1

Set the following model parameters:

temporal block resolution ["1h","1h","1h","1h","4h","4h"]
model start: 2021-05-01T00:00:00

We will then have the ramp down limit constraint with t_before and t_after as follows:

t_before=2021-05-01T03:00 (1_hour)->2021-05-01T04:00,
t_after=2021-05-01T04:00 (4_hours) ->2021-05-01T08:00

Realized constraint:

3.6 * unit_flow[t_after] >= 640 * units_on[t_after] -200 * units_on[t_before] + 0.9 * unit_flow[t_before] -800 * units_shut_down[t_after]

Let's suppose the the demand (and thus unit output) is as follows:

unit_flow[t_before] = 90
unit_flow[t_after] = 70

It should be possible within the 0.1 ramp limit but the constraint will not be fulfilled unless units_shut_down[t_after] > 0.

Expected behavior
The ramping should be possible without shutting down units.

Desktop (please complete the following information):

• Which SpineOpt version? 0.7.2
• Which Spineinterface version? 0.13.7

[nnhjy]

Looks interesting to me as well. I got the ramp_down constraint as follows:
unit_flow[t_before(1h)] - 4 unit_flow[t_after(4h)] - 200 units_on[t_before(1h)] + 640 units_on[t_after(4h)] - 800 units_shut_down[t_after(4h)] <= 0

They look similar, but why are the coefficients for the two unit_flow 1 and 4 in mine instead of 0.9 and 3.6? I didn't use ramp_up_limit = 0.1, but it should be irrelevant to the ramp_down constraint.

The results seem to illustrate some magic under the hood. The unit_shut_down line overlaps with the unit_started_up.

For the units_shut_down [t_after]> 0, the model creates a counterpart units_started_up with the same value. An interpretation of that could be "starting up and shutting down simultaneously equals no change of online status."

However, it would be better if the variable units_shut_down is automatically activated when the parameter ramp_down_limit is used.

[Mastomaki]

你 好!

I have no idea why my coefficient is 3.6 instead of 4.

But the main problem is that unit_flow[t_after] is multiplied by overlap_duration(t_after, t), in this case 4 hours, which makes it look like the unit is ramping up a lot.

Yes the units_started_up must be equal to units_shut_down because fix_units_on was set. This was also discussed in issue Simultaneous starts/shutdowns.

[Mastomaki]

My coefficient is 3.6 instead of 4 because I defined a fix_ratio_out_in_unit_flow 0.9 between fuel and output.

[Mastomaki]

Constraint code
This is my attempt to correct it:

function _build_constraint_ramp_down(m::Model, u, ng, d, s_path, t_before, t_after)
    @fetch units_on, units_shut_down, unit_flow = m.ext[:spineopt].variables
    t0 = _analysis_time(m)

    # auxiliary functions for calculating time durations
    overlap_duration_flow = t1 -> sum(overlap_duration(t1, t)
                for (u, n, d, s, t) in unit_flow_indices(m; unit=u, node=ng, 
                    direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t1));
                init=0)
    
    overlap_duration_units_on = t1 -> sum(overlap_duration(t1, t)
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t1);
                init=0)
    
    overlap_duration_switched = t1 -> sum(overlap_duration(t1, t)
                for (u, s, t) in units_switched_indices(m; unit=u, stochastic_scenario=s_path, t=t1);
                init=0)    

    @build_constraint(
        + sum(
            + unit_flow[u, n, d, s, t] * overlap_duration(t_before, t)
            for (u, n, d, s, t) in unit_flow_indices(
                m; unit=u, node=ng, direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t_before)
            )
            if !is_reserve_node(node=n);
            init=0,
        ) / overlap_duration_flow(t_before)

        - sum(
            + unit_flow[u, n, d, s, t] * overlap_duration(t_after, t)
            for (u, n, d, s, t) in unit_flow_indices(
                m; unit=u, node=ng, direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t_after)
            )
            if !is_reserve_node(node=n);
            init=0,
        ) / overlap_duration_flow(t_after)

        + sum(
            + unit_flow[u, n, d, s, t] * overlap_duration(t_after, t)
            for (u, n, d, s, t) in unit_flow_indices(
                m; unit=u, node=ng, direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t_after)
            )
            if is_reserve_node(node=n)
            && _switch(d; to_node=downward_reserve, from_node=upward_reserve)(node=n)
            && !is_non_spinning(node=n);
            init=0,
        ) / overlap_duration_flow(t_after)
        <=
         (
            + sum(  (
                    + _shut_down_limit(m, u, ng, d, s, t0, t_after)
                    - _minimum_operating_point(m, u, ng, d, s, t0, t_after)
                    - _ramp_down_limit(m, u, ng, d, s, t0, t_after)
                    )
                * _unit_flow_capacity(m, u, ng, d, s, t0, t_after)
                * units_shut_down[u, s, t]
                * duration(t)
                for (u, s, t) in units_switched_indices(m; unit=u, stochastic_scenario=s_path, t=t_after);
                init=0,
            ) / overlap_duration_switched(t_after)

            + sum(
                - _minimum_operating_point(m, u, ng, d, s, t0, t_after)
                * _unit_flow_capacity(m, u, ng, d, s, t0, t_after)
                * units_on[u, s, t]
                * duration(t)
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t_after);
                init=0,
            ) / overlap_duration_units_on(t_after)

            + sum(
                + _minimum_operating_point(m, u, ng, d, s, t0, t_after)
                * _unit_flow_capacity(m, u, ng, d, s, t0, t_after)
                * units_on[u, s, t]
                * duration(t)
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t_before);
                init=0,
            ) / overlap_duration_units_on(t_before)

            + sum(
                + _ramp_down_limit(m, u, ng, d, s, t0, t_after)
                * _unit_flow_capacity(m, u, ng, d, s, t0, t_after)
                * units_on[u, s, t]
                * duration(t)
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t_after);
                init=0,
            )
        )
        
    )
end

Results
As I result for

fix_units_on: 1
unit_capacity: 100
minimum_operating_point: 0.4
ramp_down_limit: 0.12
fix_ratio_out_in_unit_flow: 0.9

t_before=2021-05-01T03:00 (1_hour)->2021-05-01T04:00,
t_after=2021-05-01T04:00 (4_hours) ->2021-05-01T08:00

I get a very reasonable result:

0.9 * unit_flow[t_after] >= -8 * units_on[t_after] -40 * units_on[t_before] + 0.9 * unit_flow[t_before] -48 * units_shut_down[t_after]

Here the 0.9 comes just from fix_ratio_out_in_unit_flow because the output flow is automatically replaced by input flow in the constraint.

[Mastomaki]

Here is my next version.

For the example:

fix_units_on: 1
unit_capacity: 100
minimum_operating_point: 0.4
ramp_down_limit: 0.12
fix_ratio_out_in_unit_flow: 0.9

temporal block resolution ["1h","1h","1h","1h","4h","4h"]
model start: 2021-05-01T00:00:00

t_before=2021-05-01T03:00 (1_hour)->2021-05-01T04:00,
t_after=2021-05-01T04:00 (4_hours) ->2021-05-01T08:00

The result for maximum ramp down for unit flow is 30 units as expected. This comes from 6 units during t_before and 24 units during t_after. More explanations in the next message.

function _build_constraint_ramp_down(m::Model, u, ng, d, s_path, t_before, t_after)

    @fetch units_on, units_shut_down, unit_flow = m.ext[:spineopt].variables

    # auxiliary functions for calculating time durations
    overlap_duration_flow = t1 -> sum(overlap_duration(t1, t)
                for (u, n, d, s, t) in unit_flow_indices(m; unit=u, node=ng, 
                    direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t1));
                init=0)
    
    overlap_duration_units_on = t1 -> sum(overlap_duration(t1, t)
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t1);
                init=0)
    
    overlap_duration_switched = t1 -> sum(overlap_duration(t1, t)
                for (u, s, t) in units_switched_indices(m; unit=u, stochastic_scenario=s_path, t=t1);
                init=0)    

    if isdefined(Main, :Infiltrator)
        Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__)
    end
    
    @build_constraint(
        + sum(
            + unit_flow[u, n, d, s, t] * overlap_duration(t_before, t)
            for (u, n, d, s, t) in unit_flow_indices(
                m; unit=u, node=ng, direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t_before)
            )
            if !is_reserve_node(node=n);
            init=0,
        ) / overlap_duration_flow(t_before)

        - sum(
            + unit_flow[u, n, d, s, t] * overlap_duration(t_after, t)
            for (u, n, d, s, t) in unit_flow_indices(
                m; unit=u, node=ng, direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t_after)
            )
            if !is_reserve_node(node=n);
            init=0,
        ) / overlap_duration_flow(t_after)

        + sum(
            + unit_flow[u, n, d, s, t] * overlap_duration(t_after, t)
            for (u, n, d, s, t) in unit_flow_indices(
                m; unit=u, node=ng, direction=d, stochastic_scenario=s_path, t=t_overlaps_t(m; t=t_after)
            )
            if is_reserve_node(node=n)

            && _switch(d; to_node=downward_reserve, from_node=upward_reserve)(node=n)
            && !is_non_spinning(node=n);
            init=0,
        ) / overlap_duration_flow(t_after)
        <=
         (
            + sum(  (
                    + _shut_down_limit(m, u, ng, d, s, t_after)
                    - _minimum_operating_point(m, u, ng, d, s, t_after)
                    )
                * _unit_flow_capacity(m, u, ng, d, s, t_after)
                * units_shut_down[u, s, t]
                * overlap_duration(t_after, t)
                for (u, s, t) in units_switched_indices(m; unit=u, stochastic_scenario=s_path, t=t_after);
                init=0,
            ) / overlap_duration_switched(t_after)

            - sum(
                + _minimum_operating_point(m, u, ng, d, s, t_after)
                * _unit_flow_capacity(m, u, ng, d, s, t_after)
                * units_on[u, s, t]
                * overlap_duration(t_after, t)
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t_after);
                init=0,
            ) / overlap_duration_units_on(t_after)

            + sum(
                + _minimum_operating_point(m, u, ng, d, s, t_after)
                * _unit_flow_capacity(m, u, ng, d, s, t_after)
                * units_on[u, s, t]
                * overlap_duration(t_before, t)
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t_before);
                init=0,
            ) / overlap_duration_units_on(t_before)

            + sum(
                + _ramp_down_limit(m, u, ng, d, s, t_after)
                * _unit_flow_capacity(m, u, ng, d, s, t_after)
                * units_on[u, s, t]
                * overlap_duration(t_before, t)
                * 0.5
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t_before);
                init=0,
            )

            + sum(
                + _ramp_down_limit(m, u, ng, d, s, t_after)
                * _unit_flow_capacity(m, u, ng, d, s, t_after)
                * units_on[u, s, t]
                * overlap_duration(t_after, t)
                * 0.5
                for (u, s, t) in units_on_indices(m; unit=u, stochastic_scenario=s_path, t=t_after);
                init=0,
            )
        )
    )
end

[Mastomaki]

RHS terms.

Terms 1 and 2 calculate the time averages of unit_flow during t_before and t_after (v_flow,ave in the Fig 1). overlap_duration_flow() calculates the total duration of unit_flow time slices during t_before and t_after.

Fig. 1. Unit flow time average calculation.

Not sure about the third term on RHS.

LHS terms

Terms 2 and 3 on the RHS calculate the effect of units_on variable. The assumption is that the units are generating at minimum operating point level when they are shutting down (during t_before) or starting up (during t_after). In case of ramping down, the effect of units_on is positive (allows more ramping) during t_before (Fig 2). These terms also calculate the time average during t_before and t_after but as far as I know it is not really needed because units_on follows the same time slices as t_before and t_after.

Term 1 corrects the effect of shut_down_limit for plants which are shutting down during t_after. So for them it is assumed that they are generating at shut_down_limit during t_before and zero during t_after. Not sure if ramp down limit is needed in this term.

Terms 4 and 5 calculate the actual ramping potential during half of t_before and t_after. No adjustments at this point for units which are shutting down or starting up.

Fig 2. Units started up and units shut down contribute to the ramping capacity. units_started_up is not explicitly in the constraint.