2d_constraints.gms

$ontext
This file is part of Backbone.

Backbone is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

Backbone is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with Backbone.  If not, see <http://www.gnu.org/licenses/>.
$offtext


* =============================================================================
* --- Constraint Equation Definitions -----------------------------------------
* =============================================================================

* --- Energy Balance ----------------------------------------------------------

q_balance(gn(grid, node), msft(m, s, f, t)) // Energy/power balance dynamics solved using implicit Euler discretization
    ${  not p_gn(grid, node, 'boundAll')
        } ..

    // The left side of the equation is the change in the state (will be zero if the node doesn't have a state)
    + p_gn(grid, node, 'energyStoredPerUnitOfState')${gn_state(grid, node)} // Unit conversion between v_state of a particular node and energy variables (defaults to 1, but can have node based values if e.g. v_state is in Kelvins and each node has a different heat storage capacity)
        * [
            + v_state(grid, node, s, f+df_central(f,t), t)                   // The difference between current
            - v_state(grid, node, s+ds_state(grid,node,s,t), f+df(f,t+dt(t)), t+dt(t))       // ... and previous state of the node
            ]

    =E=

    // The right side of the equation contains all the changes converted to energy terms
    + p_stepLength(m, f, t) // Multiply with the length of the timestep to convert power into energy
        * (
            // Self discharge out of the model boundaries
            - p_gn(grid, node, 'selfDischargeLoss')${ gn_state(grid, node) }
                * v_state(grid, node, s, f+df_central(f,t), t) // The current state of the node

            // Energy diffusion from this node to neighbouring nodes
            - sum(gnn_state(grid, node, to_node),
                + p_gnn(grid, node, to_node, 'diffCoeff')
                    * v_state(grid, node, s, f+df_central(f,t), t)
                ) // END sum(to_node)

            // Energy diffusion from neighbouring nodes to this node
            + sum(gnn_state(grid, from_node, node),
                + p_gnn(grid, from_node, node, 'diffCoeff')
                    * v_state(grid, from_node, s, f+df_central(f,t), t) // Incoming diffusion based on the state of the neighbouring node
                ) // END sum(from_node)

            // Controlled energy transfer, applies when the current node is on the left side of the connection
            - sum(gn2n_directional(grid, node, node_),
                + (1 - p_gnn(grid, node, node_, 'transferLoss')) // Reduce transfer losses
                    * v_transfer(grid, node, node_, s, f, t)
                + p_gnn(grid, node, node_, 'transferLoss') // Add transfer losses back if transfer is from this node to another node
                    * v_transferRightward(grid, node, node_, s, f, t)
                ) // END sum(node_)

            // Controlled energy transfer, applies when the current node is on the right side of the connection
            + sum(gn2n_directional(grid, node_, node),
                + v_transfer(grid, node_, node, s, f, t)
                - p_gnn(grid, node_, node, 'transferLoss') // Reduce transfer losses if transfer is from another node to this node
                    * v_transferRightward(grid, node_, node, s, f, t)
                ) // END sum(node_)

            // Interactions between the node and its units
            + sum(gnuft(grid, node, unit, f, t),
                + v_gen(grid, node, unit, s, f, t) // Unit energy generation and consumption
                )

            // Spilling energy out of the endogenous grids in the model
            - v_spill(grid, node, s, f, t)${node_spill(node)}

            // Power inflow and outflow timeseries to/from the node
            + ts_influx_(grid, node, f, t, s)   // Incoming (positive) and outgoing (negative) absolute value time series

            // Dummy generation variables, for feasibility purposes
            + vq_gen('increase', grid, node, s, f, t) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
            - vq_gen('decrease', grid, node, s, f, t) // Note! When stateSlack is permitted, have to take caution with the penalties so that it will be used first
    ) // END * p_stepLength
;

* --- Reserve Demand ----------------------------------------------------------
// NOTE! Currently, there are multiple identical instances of the reserve balance equation being generated for each forecast branch even when the reserves are committed and identical between the forecasts.
// NOTE! This could be solved by formulating a new "ft_reserves" set to cover only the relevant forecast-time steps, but it would possibly make the reserves even more confusing.

q_resDemand(restypeDirectionNode(restype, up_down, node), sft(s, f, t))
    ${  ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')
        and not [ restypeReleasedForRealization(restype)
                  and sft_realized(s, f, t)]
        } ..

    // Reserve provision by capable units on this node
    + sum(nuft(node, unit, f, t)${nuRescapable(restype, up_down, node, unit)},
        + v_reserve(restype, up_down, node, unit, s, f+df_reserves(node, restype, f, t), t)
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision from other reserve categories when they can be shared
    + sum((nuft(node, unit, f, t), restype_)${p_nuRes2Res(node, unit, restype_, up_down, restype)},
        + v_reserve(restype_, up_down, node, unit, s, f+df_reserves(node, restype_, f, t), t)
            * p_nuRes2Res(node, unit, restype_, up_down, restype)
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                    * p_nuReserves(node, unit, restype_, 'reserveReliability')
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision to this node via transfer links
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, up_down, node_, node)},
        + (1 - p_gnn(grid, node_, node, 'transferLoss') )
            * v_resTransferRightward(restype, up_down, node_, node, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
        ) // END sum(gn2n_directional)
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, up_down, node_, node)},
        + (1 - p_gnn(grid, node, node_, 'transferLoss') )
            * v_resTransferLeftward(restype, up_down, node, node_, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
        ) // END sum(gn2n_directional)

    =G=

    // Demand for reserves
    + ts_reserveDemand_(restype, up_down, node, f, t)${p_nReserves(node, restype, 'use_time_series')}
    + p_nReserves(node, restype, up_down)${not p_nReserves(node, restype, 'use_time_series')}

    // Reserve demand increase because of units
    + sum(nuft(node, unit, f, t)${p_nuReserves(node, unit, restype, 'reserve_increase_ratio')}, // Could be better to have 'reserve_increase_ratio' separately for up and down directions
        + sum(gnu(grid, node, unit), v_gen(grid, node, unit, s, f, t)) // Reserve sets and variables are currently lacking the grid dimension...
            * p_nuReserves(node, unit, restype, 'reserve_increase_ratio')
        ) // END sum(nuft)

    // Reserve provisions to another nodes via transfer links
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, up_down, node, node_)},   // If trasferring reserves to another node, increase your own reserves by same amount
        + v_resTransferRightward(restype, up_down, node, node_, s, f+df_reserves(node, restype, f, t), t)
        ) // END sum(gn2n_directional)
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, up_down, node, node_)},   // If trasferring reserves to another node, increase your own reserves by same amount
        + v_resTransferLeftward(restype, up_down, node_, node, s, f+df_reserves(node, restype, f, t), t)
        ) // END sum(gn2n_directional)

    // Reserve demand feasibility dummy variables
    - vq_resDemand(restype, up_down, node, s, f+df_reserves(node, restype, f, t), t)
    - vq_resMissing(restype, up_down, node, s, f+df_reserves(node, restype, f, t), t)${ft_reservesFixed(node, restype, f+df_reserves(node, restype, f, t), t)}
;

* --- N-1 Reserve Demand ----------------------------------------------------------
// NOTE! Currently, there are multiple identical instances of the reserve balance equation being generated for each forecast branch even when the reserves are committed and identical between the forecasts.
// NOTE! This could be solved by formulating a new "ft_reserves" set to cover only the relevant forecast-time steps, but it would possibly make the reserves even more confusing.

q_resDemandLargestInfeedUnit(grid, restypeDirectionNode(restype, 'up', node), unit_fail(unit_), sft(s, f, t))
    ${  ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')
        and gn(grid, node)
        and not [ restypeReleasedForRealization(restype)
            and ft_realized(f, t)
            ]
        and p_nuReserves(node, unit_, restype, 'portion_of_infeed_to_reserve')
        } ..

    // Reserve provision by capable units on this node excluding the failing one
    + sum(nuft(node, unit, f, t)${nuRescapable(restype, 'up', node, unit) and (ord(unit_) ne ord(unit))},
        + v_reserve(restype, 'up', node, unit, s, f+df_reserves(node, restype, f, t), t)
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision from other reserve categories when they can be shared
    + sum((nuft(node, unit, f, t), restype_)${p_nuRes2Res(node, unit, restype_, 'up', restype)},
        + v_reserve(restype_, 'up', node, unit, s, f+df_reserves(node, restype_, f, t), t)
            * p_nuRes2Res(node, unit, restype_, 'up', restype)
            * [ // Account for reliability of reserves
                + 1${sft_realized(s, f+df_reserves(node, restype, f, t), t)} // reserveReliability limits the reliability of reserves locked ahead of time.
                + p_nuReserves(node, unit, restype, 'reserveReliability')${not sft_realized(s, f+df_reserves(node, restype, f, t), t)}
                    * p_nuReserves(node, unit, restype_, 'reserveReliability')
                ] // END * v_reserve
        ) // END sum(nuft)

    // Reserve provision to this node via transfer links
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'up', node_, node)},
        + (1 - p_gnn(grid, node_, node, 'transferLoss') )
            * v_resTransferRightward(restype, 'up', node_, node, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
        ) // END sum(gn2n_directional)
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'up', node_, node)},
        + (1 - p_gnn(grid, node, node_, 'transferLoss') )
            * v_resTransferLeftward(restype, 'up', node, node_, s, f+df_reserves(node_, restype, f, t), t) // Reserves from another node - reduces the need for reserves in the node
        ) // END sum(gn2n_directional)

    =G=

    // Demand for reserves of the failing one
    v_gen(grid,node,unit_,s,f,t) * p_nuReserves(node, unit_, restype, 'portion_of_infeed_to_reserve')

    // Reserve provisions to another nodes via transfer links
    + sum(gn2n_directional(grid, node, node_)${restypeDirectionNodeNode(restype, 'up', node, node_)},   // If trasferring reserves to another node, increase your own reserves by same amount
        + v_resTransferRightward(restype, 'up', node, node_, s, f+df_reserves(node, restype, f, t), t)
        ) // END sum(gn2n_directional)
    + sum(gn2n_directional(grid, node_, node)${restypeDirectionNodeNode(restype, 'up', node, node_)},   // If trasferring reserves to another node, increase your own reserves by same amount
        + v_resTransferLeftward(restype, 'up', node_, node, s, f+df_reserves(node, restype, f, t), t)
        ) // END sum(gn2n_directional)

    // Reserve demand feasibility dummy variables
    - vq_resDemand(restype, 'up', node, s, f+df_reserves(node, restype, f, t), t)
    - vq_resMissing(restype, 'up', node, s, f+df_reserves(node, restype, f, t), t)${ft_reservesFixed(node, restype, f+df_reserves(node, restype, f, t), t)}
;
* --- Maximum Downward Capacity -----------------------------------------------

q_maxDownward(gnu(grid, node, unit), msft(m, s, f, t))
    ${  gnuft(grid, node, unit, f, t)
        and {
            [   ord(t) < tSolveFirst + smax(restype, p_nReserves(node, restype, 'reserve_length')) // Unit is either providing
                and sum(restype, nuRescapable(restype, 'down', node, unit)) // downward reserves
                ]
            // NOTE!!! Could be better to form a gnuft_reserves subset?
            or [ // the unit has an online variable
                uft_online(unit, f, t)
                and [
                    (unit_minLoad(unit) and p_gnu(grid, node, unit, 'unitSizeGen')) // generators with a min. load
                    or p_gnu(grid, node, unit, 'maxCons') // or consuming units with an online variable
                    ]
                ] // END or
            or [ // consuming units with investment possibility
                gnu_input(grid, node, unit)
                and [unit_investLP(unit) or unit_investMIP(unit)]
                ]
        }} ..

    // Energy generation/consumption
    + v_gen(grid, node, unit, s, f, t)

    // Considering output constraints (e.g. cV line)
    + sum(gngnu_constrainedOutputRatio(grid, node, grid_output, node_, unit),
        + p_gnu(grid_output, node_, unit, 'cV')
            * v_gen(grid_output, node_, unit, s, f, t)
        ) // END sum(gngnu_constrainedOutputRatio)

    // Downward reserve participation
    - sum(nuRescapable(restype, 'down', node, unit)${ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')},
        + v_reserve(restype, 'down', node, unit, s, f+df_reserves(node, restype, f, t), t) // (v_reserve can be used only if the unit is capable of providing a particular reserve)
        ) // END sum(nuRescapable)

    =G= // Must be greater than minimum load or maximum consumption  (units with min-load and both generation and consumption are not allowed)

    // Generation units, greater than minload
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
            + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
            + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
            ) // END sum(effGroup)
        * [ // Online variables should only be generated for units with restrictions
            + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f+df_central(f,t), t)} // LP online variant
            + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f+df_central(f,t), t)} // MIP online variant
            ] // END v_online

    // Units in run-up phase neet to keep up with the run-up rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(unitStarttype(unit, starttype)${uft_startupTrajectory(unit, f, t)},
            sum(runUpCounter(unit, counter)${t_active(t+dt_trajectory(counter))}, // Sum over the run-up intervals
                + [
                    + v_startup_LP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    + v_startup_MIP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    ]
                    * p_uCounter_runUpMin(unit, counter)
                ) // END sum(runUpCounter)
            ) // END sum(unitStarttype)

    // Units in shutdown phase need to keep up with the shutdown rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(shutdownCounter(unit, counter)${t_active(t+dt_trajectory(counter)) and uft_shutdownTrajectory(unit, f, t)}, // Sum over the shutdown intervals
            + [
                + v_shutdown_LP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                + v_shutdown_MIP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                ]
                * p_uCounter_shutdownMin(unit, counter)
            ) // END sum(shutdownCounter)

    // Consuming units, greater than maxCons
    // Available capacity restrictions
    - p_unit(unit, 'availability')
        * [
            // Capacity factors for flow units
            + sum(flowUnit(flow, unit),
                + ts_cf_(flow, node, f, t, s)
                ) // END sum(flow)
            + 1${not unit_flow(unit)}
            ] // END * p_unit(availability)
        * [
            // Online capacity restriction
            + p_gnu(grid, node, unit, 'maxCons')${not uft_online(unit, f, t)} // Use initial maximum if no online variables
            + p_gnu(grid, node, unit, 'unitSizeCons')
                * [
                    // Capacity online
                    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
                    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}

                    // Investments to additional non-online capacity
                    + sum(t_invest(t_)${    ord(t_)<=ord(t)
                                            and not uft_online(unit, f, t)
                                            },
                        + v_invest_LP(unit, t_)${unit_investLP(unit)} // NOTE! v_invest_LP also for consuming units is positive
                        + v_invest_MIP(unit, t_)${unit_investMIP(unit)} // NOTE! v_invest_MIP also for consuming units is positive
                        ) // END sum(t_invest)
                    ] // END * p_gnu(unitSizeCons)
            ] // END * p_unit(availability)
;

* --- Maximum Upwards Capacity ------------------------------------------------

q_maxUpward(gnu(grid, node, unit), msft(m, s, f, t))
    ${  gnuft(grid, node, unit, f, t)
        and {
            [   ord(t) < tSolveFirst + smax(restype, p_nReserves(node, restype, 'reserve_length')) // Unit is either providing
                and sum(restype, nuRescapable(restype, 'up', node, unit)) // upward reserves
                ]
            or [
                uft_online(unit, f, t) // or the unit has an online variable
                and [
                    [unit_minLoad(unit) and p_gnu(grid, node, unit, 'unitSizeCons')] // consuming units with min_load
                    or [p_gnu(grid, node, unit, 'maxGen')]                          // generators with an online variable
                    ]
                ]
            or [
                gnu_output(grid, node, unit) // generators with investment possibility
                and (unit_investLP(unit) or unit_investMIP(unit))
                ]
        }}..

    // Energy generation/consumption
    + v_gen(grid, node, unit, s, f, t)

    // Considering output constraints (e.g. cV line)
    + sum(gngnu_constrainedOutputRatio(grid, node, grid_output, node_, unit),
        + p_gnu(grid_output, node_, unit, 'cV')
            * v_gen(grid_output, node_, unit, s, f, t)
        ) // END sum(gngnu_constrainedOutputRatio)

    // Upwards reserve participation
    + sum(nuRescapable(restype, 'up', node, unit)${ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')},
        + v_reserve(restype, 'up', node, unit, s, f+df_reserves(node, restype, f, t), t)
        ) // END sum(nuRescapable)

    =L= // must be less than available/online capacity

    // Consuming units
    - p_gnu(grid, node, unit, 'unitSizeCons')
        * sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
            + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
            + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
            ) // END sum(effGroup)
        * [
            + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)} // Consuming units are restricted by their min. load (consuming is negative)
            + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)} // Consuming units are restricted by their min. load (consuming is negative)
            ] // END * p_gnu(unitSizeCons)

    // Generation units
    // Available capacity restrictions
    + p_unit(unit, 'availability') // Generation units are restricted by their (available) capacity
        * [
            // Capacity factor for flow units
            + sum(flowUnit(flow, unit),
                + ts_cf_(flow, node, f, t, s)
                ) // END sum(flow)
            + 1${not unit_flow(unit)}
            ] // END * p_unit(availability)
        * [
            // Online capacity restriction
            + p_gnu(grid, node, unit, 'maxGen')${not uft_online(unit, f, t)} // Use initial maxGen if no online variables
            + p_gnu(grid, node, unit, 'unitSizeGen')
                * [
                    // Capacity online
                    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f ,t)}
                    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}

                    // Investments to non-online capacity
                    + sum(t_invest(t_)${    ord(t_)<=ord(t)
                                            and not uft_online(unit, f ,t)
                                            },
                        + v_invest_LP(unit, t_)${unit_investLP(unit)}
                        + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
                        ) // END sum(t_invest)
                    ] // END * p_gnu(unitSizeGen)
            ] // END * p_unit(availability)

    // Units in run-up phase neet to keep up with the run-up rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(unitStarttype(unit, starttype)${uft_startupTrajectory(unit, f, t)},
            sum(runUpCounter(unit, counter)${t_active(t+dt_trajectory(counter))}, // Sum over the run-up intervals
                + [
                    + v_startup_LP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    + v_startup_MIP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    ]
                    * p_uCounter_runUpMax(unit, counter)
                ) // END sum(runUpCounter)
            ) // END sum(unitStarttype)

    // Units in shutdown phase need to keep up with the shutdown rate
    + p_gnu(grid, node, unit, 'unitSizeGen')
        * sum(shutdownCounter(unit, counter)${t_active(t+dt_trajectory(counter)) and uft_shutdownTrajectory(unit, f, t)}, // Sum over the shutdown intervals
            + [
                + v_shutdown_LP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                + v_shutdown_MIP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                    ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                ]
                * p_uCounter_shutdownMax(unit, counter)
            ) // END sum(shutdownCounter)
;

* --- Reserve Provision of Units with Investments -----------------------------

q_reserveProvision(nuRescapable(restypeDirectionNode(restype, up_down, node), unit), sft(s, f, t))
    ${  ord(t) <= tSolveFirst + p_nReserves(node, restype, 'reserve_length')
        and nuft(node, unit, f, t)
        and (unit_investLP(unit) or unit_investMIP(unit))
        and not ft_reservesFixed(node, restype, f+df_reserves(node, restype, f, t), t)
        } ..

    + v_reserve(restype, up_down, node, unit, s, f+df_reserves(node, restype, f, t), t)

    =L=

    + p_nuReserves(node, unit, restype, up_down)
        * [
            + sum(grid, p_gnu(grid, node, unit, 'maxGen') + p_gnu(grid, node, unit, 'maxCons') )  // Reserve sets and variables are currently lacking the grid dimension...
            + sum(t_invest(t_)${ ord(t_)<=ord(t) },
                + v_invest_LP(unit, t_)${unit_investLP(unit)}
                    * sum(grid, p_gnu(grid, node, unit, 'unitSizeTot')) // Reserve sets and variables are currently lacking the grid dimension...
                + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
                    * sum(grid, p_gnu(grid, node, unit, 'unitSizeTot')) // Reserve sets and variables are currently lacking the grid dimension...
                ) // END sum(t_)
            ]
        * p_unit(unit, 'availability') // Taking into account availability...
        * [
            // ... and capacity factor for flow units
            + sum(flowUnit(flow, unit),
                + ts_cf_(flow, node, f, t, s)
                ) // END sum(flow)
            + 1${not unit_flow(unit)}
            ] // How to consider reserveReliability in the case of investments when we typically only have "realized" time steps?
;

* --- Unit Startup and Shutdown -----------------------------------------------

q_startshut(ms(m, s), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        }..

    // Units currently online
    + v_online_LP (unit, s, f+df_central(f,t), t)${uft_onlineLP (unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}

    // Units previously online
    // The same units
    - v_online_LP (unit, s+ds(s,t), f+df(f,t+dt(t)), t+dt(t))${ uft_onlineLP_withPrevious(unit, f+df(f,t+dt(t)), t+dt(t))
                                                             and not uft_aggregator_first(unit, f, t) } // This reaches to tFirstSolve when dt = -1
    - v_online_MIP(unit, s+ds(s,t), f+df(f,t+dt(t)), t+dt(t))${ uft_onlineMIP_withPrevious(unit, f+df(f,t+dt(t)), t+dt(t))
                                                             and not uft_aggregator_first(unit, f, t) }

    // Aggregated units just before they are turned into aggregator units
    - sum(unit_${unitAggregator_unit(unit, unit_)},
        + v_online_LP (unit_, s, f+df(f,t+dt(t)), t+dt(t))${uft_onlineLP_withPrevious(unit_, f+df(f,t+dt(t)), t+dt(t))}
        + v_online_MIP(unit_, s, f+df(f,t+dt(t)), t+dt(t))${uft_onlineMIP_withPrevious(unit_, f+df(f,t+dt(t)), t+dt(t))}
        )${uft_aggregator_first(unit, f, t)} // END sum(unit_)

    =E=

    // Unit startup and shutdown

    // Add startup of units dt_toStartup before the current t (no start-ups for aggregator units before they become active)
    + sum(unitStarttype(unit, starttype),
        + v_startup_LP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t))
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
        + v_startup_MIP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t))
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
        )${not [unit_aggregator(unit) and ord(t) + dt_toStartup(unit, t) <= tSolveFirst + p_unit(unit, 'lastStepNotAggregated')]} // END sum(starttype)

    // NOTE! According to 3d_setVariableLimits,
    // cannot start a unit if the time when the unit would become online is outside
    // the horizon when the unit has an online variable
    // --> no need to add start-ups of aggregated units to aggregator units

    // Shutdown of units at time t
    - v_shutdown_LP(unit, s, f, t)
        ${ uft_onlineLP(unit, f, t) }
    - v_shutdown_MIP(unit, s, f, t)
        ${ uft_onlineMIP(unit, f, t) }
;

*--- Startup Type -------------------------------------------------------------
// !!! NOTE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// This formulation doesn't work as intended when unitCount > 1, as one recent
// shutdown allows for multiple hot/warm startups on subsequent time steps.
// Pending changes.

q_startuptype(ms(m, s), starttypeConstrained(starttype), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        and unitStarttype(unit, starttype)
        } ..

    // Startup type
    + v_startup_LP(unit, starttype, s, f, t)${ uft_onlineLP(unit, f, t) }
    + v_startup_MIP(unit, starttype, s, f, t)${ uft_onlineMIP(unit, f, t) }

    =L=

    // Subunit shutdowns within special startup timeframe
    + sum(unitCounter(unit, counter)${  dt_starttypeUnitCounter(starttype, unit, counter)
                                        and t_active(t+(dt_starttypeUnitCounter(starttype, unit, counter)+1))
                                        },
        + v_shutdown_LP(unit, s, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1))
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)) }
        + v_shutdown_MIP(unit, s, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1))
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)), t+(dt_starttypeUnitCounter(starttype, unit, counter)+1)) }
        ) // END sum(counter)

    // NOTE: for aggregator units, shutdowns for aggregated units are not considered
;


*--- Online Limits with Startup Type Constraints and Investments --------------

q_onlineLimit(ms(m, s), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        and {
            p_unit(unit, 'minShutdownHours')
            or p_u_runUpTimeIntervals(unit)
            or unit_investLP(unit)
            or unit_investMIP(unit)
        }} ..

    // Online variables
    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f ,t)}

    =L=

    // Number of existing units
    + p_unit(unit, 'unitCount')

    // Number of units unable to become online due to restrictions
    - sum(unitCounter(unit, counter)${  dt_downtimeUnitCounter(unit, counter)
                                        and t_active(t+(dt_downtimeUnitCounter(unit, counter) + 1))
                                        },
        + v_shutdown_LP(unit, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
        + v_shutdown_MIP(unit, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
        ) // END sum(counter)

    // Number of units unable to become online due to restrictions (aggregated units in the past horizon or if they have an online variable)
    - sum(unitAggregator_unit(unit, unit_),
        + sum(unitCounter(unit, counter)${  dt_downtimeUnitCounter(unit, counter)
                                            and t_active(t+(dt_downtimeUnitCounter(unit, counter) + 1))
                                            },
            + v_shutdown_LP(unit_, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
                ${ uft_onlineLP_withPrevious(unit_, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
            + v_shutdown_MIP(unit_, s, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1))
                ${ uft_onlineMIP_withPrevious(unit_, f+df(f,t+(dt_downtimeUnitCounter(unit, counter) + 1)), t+(dt_downtimeUnitCounter(unit, counter) + 1)) }
            ) // END sum(counter)
        )${unit_aggregator(unit)} // END sum(unit_)

    // Investments into units
    + sum(t_invest(t_)${ord(t_)<=ord(t)},
        + v_invest_LP(unit, t_)${unit_investLP(unit)}
        + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
        ) // END sum(t_invest)
;

*--- Both q_offlineAfterShutdown and q_onlineOnStartup work when there is only one unit.
*    These equations prohibit single units turning on and off at the same time step.
*    Unfortunately there seems to be no way to prohibit this when unit count is > 1.
*    (it shouldn't be worthwhile anyway if there is a startup cost, but it can fall within the solution gap).
q_onlineOnStartUp(s_active(s), uft_online(unit, f, t))
    ${  sft(s, f, t)
        and sum(starttype, unitStarttype(unit, starttype))
        }..

    // Units currently online
    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}

    =G=

    + sum(unitStarttype(unit, starttype),
        + v_startup_LP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) //dt_toStartup displaces the time step to the one where the unit would be started up in order to reach online at t
            ${ uft_onlineLP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
        + v_startup_MIP(unit, starttype, s, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) //dt_toStartup displaces the time step to the one where the unit would be started up in order to reach online at t
            ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+dt_toStartup(unit, t)), t+dt_toStartup(unit, t)) }
      ) // END sum(starttype)
;

q_offlineAfterShutdown(s_active(s), uft_online(unit, f, t))
    ${  sft(s, f, t)
        and sum(starttype, unitStarttype(unit, starttype))
        }..

    // Number of existing units
    + p_unit(unit, 'unitCount')

    // Investments into units
    + sum(t_invest(t_)${ord(t_)<=ord(t)},
        + v_invest_LP(unit, t_)${unit_investLP(unit)}
        + v_invest_MIP(unit, t_)${unit_investMIP(unit)}
        ) // END sum(t_invest)

    // Units currently online
    - v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    - v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}

    =G=

    + v_shutdown_LP(unit, s, f, t)
        ${ uft_onlineLP(unit, f, t) }
    + v_shutdown_MIP(unit, s, f, t)
        ${ uft_onlineMIP(unit, f, t) }
;

*--- Minimum Unit Uptime ------------------------------------------------------

q_onlineMinUptime(ms(m, s), uft_online(unit, f, t))
    ${  msft(m, s, f, t)
        and  p_unit(unit, 'minOperationHours')
        } ..

    // Units currently online
    + v_online_LP(unit, s, f+df_central(f,t), t)${uft_onlineLP(unit, f, t)}
    + v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}

    =G=

    // Units that have minimum operation time requirements active
    + sum(unitCounter(unit, counter)${  dt_uptimeUnitCounter(unit, counter)
                                        and t_active(t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) // Don't sum over counters that don't point to an active time step
                                        },
        + sum(unitStarttype(unit, starttype),
            + v_startup_LP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
                ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
            + v_startup_MIP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
                ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
            ) // END sum(starttype)
        ) // END sum(counter)

    // Units that have minimum operation time requirements active (aggregated units in the past horizon or if they have an online variable)
    + sum(unitAggregator_unit(unit, unit_),
        + sum(unitCounter(unit, counter)${  dt_uptimeUnitCounter(unit, counter)
                                            and t_active(t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) // Don't sum over counters that don't point to an active time step
                                            },
            + sum(unitStarttype(unit, starttype),
                + v_startup_LP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
                    ${ uft_onlineLP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
                + v_startup_MIP(unit, starttype, s, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1))
                    ${ uft_onlineMIP_withPrevious(unit, f+df(f,t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)), t+(dt_uptimeUnitCounter(unit, counter)+dt_toStartup(unit, t) + 1)) }
                ) // END sum(starttype)
            ) // END sum(counter)
        )${unit_aggregator(unit)} // END sum(unit_)
;

* --- Cyclic Boundary Conditions for Online State -----------------------------

q_onlineCyclic(uss_bound(unit, s_, s), m)
    ${  ms(m, s_)
        and ms(m, s)
        and tSolveFirst = mSettings(m, 't_start')
        }..

    // Initial value of the state of the unit at the start of the sample
    + sum(mst_start(m, s, t),
        + sum(sft(s, f, t),
            + v_online_LP(unit, s, f+df(f,t+dt(t)), t+dt(t))
                ${uft_onlineLP_withPrevious(unit, f+df(f,t+dt(t)), t+dt(t))}
            + v_online_MIP(unit, s, f+df(f,t+dt(t)), t+dt(t))
                ${uft_onlineMIP_withPrevious(unit, f+df(f,t+dt(t)), t+dt(t))}
            ) // END sum(ft)
        ) // END sum(mst_start)

    =E=

    // State of the unit at the end of the sample
    + sum(mst_end(m, s_, t_),
        + sum(sft(s_, f_, t_),
            + v_online_LP(unit, s_, f_, t_)${uft_onlineLP(unit, f_, t_)}
            + v_online_MIP(unit, s_, f_, t_)${uft_onlineMIP(unit, f_, t_)}
            ) // END sum(ft)
        ) // END sum(mst_end)
;

* --- Ramp Constraints --------------------------------------------------------

q_genRamp(ms(m, s), gnuft_ramp(grid, node, unit, f, t))
    ${  ord(t) > msStart(m, s) + 1
        and msft(m, s, f, t)
        } ..

    + v_genRamp(grid, node, unit, s, f, t)
        * p_stepLength(m, f, t)

    =E=

    // Change in generation over the interval: v_gen(t) - v_gen(t-1)
    + v_gen(grid, node, unit, s, f, t)

    // Unit generation at t-1 (except aggregator units right before the aggregation threshold, see next term)
    - v_gen(grid, node, unit, s+ds(s,t), f+df(f,t+dt(t)), t+dt(t))${not uft_aggregator_first(unit, f, t)}
    // Unit generation at t-1, aggregator units right before the aggregation threshold
    + sum(unit_${unitAggregator_unit(unit, unit_)},
        - v_gen(grid, node, unit_, s+ds(s,t), f+df(f,t+dt(t)), t+dt(t))
      )${uft_aggregator_first(unit, f, t)}
;

* --- Ramp Up Limits ----------------------------------------------------------

q_rampUpLimit(ms(m, s), gnuft_ramp(grid, node, unit, f, t))
    ${  ord(t) > msStart(m, s) + 1
        and msft(m, s, f, t)
        and p_gnu(grid, node, unit, 'maxRampUp')
        and [ sum(restype, nuRescapable(restype, 'up', node, unit))
              or uft_online(unit, f, t)
              or unit_investLP(unit)
              or unit_investMIP(unit)
              ]
        } ..

    // Ramp speed of the unit?
    + v_genRamp(grid, node, unit, s, f, t)
    + sum(nuRescapable(restype, 'up', node, unit)${ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')},
        + v_reserve(restype, 'up', node, unit, s, f+df_reserves(node, restype, f, t), t) // (v_reserve can be used only if the unit is capable of providing a particular reserve)
        ) // END sum(nuRescapable)
        / p_stepLength(m, f, t)

    =L=

    // Ramping capability of units without an online variable
    + (
        + ( p_gnu(grid, node, unit, 'maxGen') + p_gnu(grid, node, unit, 'maxCons') )${not uft_online(unit, f, t)}
        + sum(t_invest(t_)${ ord(t_)<=ord(t) },
            + v_invest_LP(unit, t_)${not uft_onlineLP(unit, f, t) and unit_investLP(unit)}
                * p_gnu(grid, node, unit, 'unitSizeTot')
            + v_invest_MIP(unit, t_)${not uft_onlineMIP(unit, f, t) and unit_investMIP(unit)}
                * p_gnu(grid, node, unit, 'unitSizeTot')
          )
      )
        * p_gnu(grid, node, unit, 'maxRampUp')
        * 60   // Unit conversion from [p.u./min] to [p.u./h]

    // Ramping capability of units with an online variable
    + (
        + v_online_LP(unit, s, f+df_central(f,t), t)
            ${uft_onlineLP(unit, f, t)}
        + v_online_MIP(unit, s, f+df_central(f,t), t)
            ${uft_onlineMIP(unit, f, t)}
      )
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * p_gnu(grid, node, unit, 'maxRampUp')
        * 60   // Unit conversion from [p.u./min] to [p.u./h]

    // Generation units not be able to ramp from zero to min. load within one time interval according to their maxRampUp
    + sum(unitStarttype(unit, starttype)${   uft_online(unit, f, t)
                                             and gnu_output(grid, node, unit)
                                             and not uft_startupTrajectory(unit, f, t)
                                             and ( + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
                                                       + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
                                                       + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
                                                     ) // END sum(effGroup)
                                                       / p_stepLength(m, f, t)
                                                   - p_gnu(grid, node, unit, 'maxRampUp')
                                                       * 60 > 0
                                                   )
                                             },
        + v_startup_LP(unit, starttype, s, f, t)
            ${ uft_onlineLP(unit, f, t) }
        + v_startup_MIP(unit, starttype, s, f, t)
            ${ uft_onlineMIP(unit, f, t) }
      ) // END sum(starttype)
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * (
            + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
                + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
                + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
              ) // END sum(effGroup)
                / p_stepLength(m, f, t)
            - p_gnu(grid, node, unit, 'maxRampUp')
                * 60   // Unit conversion from [p.u./min] to [p.u./h]
          ) // END * v_startup

    // Units in the run-up phase need to keep up with the run-up rate
    + p_gnu(grid, node, unit, 'unitSizeTot')
        * sum(unitStarttype(unit, starttype)${uft_startupTrajectory(unit, f, t)},
            sum(runUpCounter(unit, counter)${t_active(t+dt_trajectory(counter))}, // Sum over the run-up intervals
                + [
                    + v_startup_LP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    + v_startup_MIP(unit, starttype, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    ]
                    * [
                        + p_unit(unit, 'rampSpeedToMinLoad')
                        + ( p_gnu(grid, node, unit, 'maxRampUp') - p_unit(unit, 'rampSpeedToMinLoad') )${ not runUpCounter(unit, counter+1) } // Ramp speed adjusted for the last run-up interval
                            * ( p_u_runUpTimeIntervalsCeil(unit) - p_u_runUpTimeIntervals(unit) )
                        ]
                    * 60 // Unit conversion from [p.u./min] into [p.u./h]
                ) // END sum(runUpCounter)
            ) // END sum(unitStarttype)

    // Shutdown of consumption units according to maxRampUp
    + [
        + v_shutdown_LP(unit, s, f, t)
            ${uft_onlineLP(unit, f, t) and gnu_input(grid, node, unit)}
        + v_shutdown_MIP(unit, s, f, t)
            ${uft_onlineMIP(unit, f, t) and gnu_input(grid, node, unit)}
        ]
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * p_gnu(grid, node, unit, 'maxRampUp')
        * 60   // Unit conversion from [p.u./min] to [p.u./h]
    // Consumption units not be able to ramp from min. load to zero within one time interval according to their maxRampUp
    + [
        + v_shutdown_LP(unit, s, f, t)
            ${ uft_onlineLP(unit, f, t) }
        + v_shutdown_MIP(unit, s, f, t)
            ${ uft_onlineMIP(unit, f, t) }
        ]
        ${  gnu_input(grid, node, unit)
            and ( + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
                      + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
                      + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
                      ) // END sum(effGroup)
                      / p_stepLength(m, f, t)
                  - p_gnu(grid, node, unit, 'maxRampUp')
                      * 60 > 0
                  )
            }
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * (
            + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
                + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
                + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
                ) // END sum(effGroup)
                / p_stepLength(m, f, t)
            - p_gnu(grid, node, unit, 'maxRampUp')
                * 60   // Unit conversion from [p.u./min] to [p.u./h]
          ) // END * v_shutdown
;

* --- Ramp Down Limits --------------------------------------------------------

q_rampDownLimit(ms(m, s), gnuft_ramp(grid, node, unit, f, t))
    ${  ord(t) > msStart(m, s) + 1
        and msft(m, s, f, t)
        and p_gnu(grid, node, unit, 'maxRampDown')
        and [ sum(restype, nuRescapable(restype, 'down', node, unit))
              or uft_online(unit, f, t)
              or unit_investLP(unit)
              or unit_investMIP(unit)
              ]
        } ..

    // Ramp speed of the unit?
    + v_genRamp(grid, node, unit, s, f, t)
    - sum(nuRescapable(restype, 'down', node, unit)${ord(t) < tSolveFirst + p_nReserves(node, restype, 'reserve_length')},
        + v_reserve(restype, 'down', node, unit, s, f+df_reserves(node, restype, f, t), t) // (v_reserve can be used only if the unit is capable of providing a particular reserve)
        ) // END sum(nuRescapable)
        / p_stepLength(m, f, t)

    =G=

    // Ramping capability of units without online variable
    - (
        + ( p_gnu(grid, node, unit, 'maxGen') + p_gnu(grid, node, unit, 'maxCons') )
            ${not uft_online(unit, f, t)}
        + sum(t_invest(t_)${ ord(t_)<=ord(t) },
            + v_invest_LP(unit, t_)
                ${not uft_onlineLP(unit, f, t) and unit_investLP(unit)}
                * p_gnu(grid, node, unit, 'unitSizeTot')
            + v_invest_MIP(unit, t_)
                ${not uft_onlineMIP(unit, f, t) and unit_investMIP(unit)}
                * p_gnu(grid, node, unit, 'unitSizeTot')
          )
      )
        * p_gnu(grid, node, unit, 'maxRampDown')
        * 60   // Unit conversion from [p.u./min] to [p.u./h]

    // Ramping capability of units that are online
    - (
        + v_online_LP(unit, s, f+df_central(f,t), t)
            ${uft_onlineLP(unit, f, t)}
        + v_online_MIP(unit, s, f+df_central(f,t), t)
            ${uft_onlineMIP(unit, f, t)}
      )
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * p_gnu(grid, node, unit, 'maxRampDown')
        * 60   // Unit conversion from [p.u./min] to [p.u./h]

    // Shutdown of generation units according to maxRampDown
    - [
        + v_shutdown_LP(unit, s, f, t)
            ${  uft_onlineLP(unit, f, t) }
        + v_shutdown_MIP(unit, s, f, t)
            ${  uft_onlineMIP(unit, f, t) }
        ]
        ${  gnu_output(grid, node, unit)
            and not uft_shutdownTrajectory(unit, f, t)
            }
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * p_gnu(grid, node, unit, 'maxRampDown')
        * 60   // Unit conversion from [p.u./min] to [p.u./h]
    // Generation units not be able to ramp from min. load to zero within one time interval according to their maxRampDown
    - [
        + v_shutdown_LP(unit, s, f, t)
            ${  uft_onlineLP(unit, f, t) }
        + v_shutdown_MIP(unit, s, f, t)
            ${  uft_onlineMIP(unit, f, t) }
        ]
        ${  gnu_output(grid, node, unit)
            and not uft_shutdownTrajectory(unit, f, t)
            and ( + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
                      + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
                      + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
                    ) // END sum(effGroup)
                    / p_stepLength(m, f, t)
                  - p_gnu(grid, node, unit, 'maxRampDown')
                      * 60 > 0
                )
        }
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * (
            + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
                + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
                + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
                ) // END sum(effGroup)
                / p_stepLength(m, f, t)
            - p_gnu(grid, node, unit, 'maxRampDown')
                * 60   // Unit conversion from [p.u./min] to [p.u./h]
          ) // END * v_shutdown

    // Units in shutdown phase need to keep up with the shutdown ramp rate
    - p_gnu(grid, node, unit, 'unitSizeGen')
        * [
            + sum(shutdownCounter(unit, counter)${t_active(t+dt_trajectory(counter)) and uft_shutdownTrajectory(unit, f, t)}, // Sum over the shutdown intervals
                + [
                    + v_shutdown_LP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    + v_shutdown_MIP(unit, s, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter))
                        ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_trajectory(counter)), t+dt_trajectory(counter)) }
                    ]
                    * [
                        + p_gnu(grid, node, unit, 'maxRampDown')${ not shutdownCounter(unit, counter-1) } // Normal maxRampDown limit applies to the time interval when v_shutdown happens, i.e. over the change from online to offline (symmetrical to v_startup)
                        + p_unit(unit, 'rampSpeedFromMinLoad')${ shutdownCounter(unit, counter-1) } // Normal trajectory ramping
                        + ( p_gnu(grid, node, unit, 'maxRampDown') - p_unit(unit, 'rampSpeedFromMinLoad') )${ shutdownCounter(unit, counter-1) and not shutdownCounter(unit, counter-2) } // Ramp speed adjusted for the first shutdown interval
                            * ( p_u_shutdownTimeIntervalsCeil(unit) - p_u_shutdownTimeIntervals(unit) )
                        ]
                ) // END sum(shutdownCounter)
            // Units need to be able to shut down after shut down trajectory
            + [
                + v_shutdown_LP(unit, s, f+df(f, t+dt_toShutdown(unit, t)), t+dt_toShutdown(unit, t))
                    ${ uft_onlineLP_withPrevious(unit, f+df(f, t+dt_toShutdown(unit, t)), t+dt_toShutdown(unit, t)) }
                + v_shutdown_MIP(unit, s, f+df(f, t+dt_toShutdown(unit, t)), t+dt_toShutdown(unit, t))
                    ${ uft_onlineMIP_withPrevious(unit, f+df(f, t+dt_toShutdown(unit, t)), t+dt_toShutdown(unit, t)) }
                ]
                ${uft_shutdownTrajectory(unit, f, t)}
                * [
                    + p_unit(unit, 'rampSpeedFromMinload')
                    + ( p_gnu(grid, node, unit, 'maxRampDown') - p_unit(unit, 'rampSpeedFromMinLoad') )${ sum(shutdownCounter(unit, counter), 1) = 1 } // Ramp speed adjusted if the unit has only one shutdown interval
                        * ( p_u_shutdownTimeIntervalsCeil(unit) - p_u_shutdownTimeIntervals(unit) )
                    ]
            ]
        * 60 // Unit conversion from [p.u./min] to [p.u./h]

    // Consumption units not be able to ramp from zero to min. load within one time interval according to their maxRampDown
    - sum(unitStarttype(unit, starttype)${   uft_online(unit, f, t)
                                             and gnu_input(grid, node, unit)
                                             and ( + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation
                                                       + p_effGroupUnit(effGroup, unit, 'lb')${not ts_effGroupUnit(effGroup, unit, 'lb', f, t)}
                                                       + ts_effGroupUnit(effGroup, unit, 'lb', f, t)
                                                     ) // END sum(effGroup)
                                                       / p_stepLength(m, f, t)
                                                   - p_gnu(grid, node, unit, 'maxRampDown')
                                                       * 60 > 0
                                                   )
                                             },
        + v_startup_LP(unit, starttype, s, f, t)
            ${ uft_onlineLP(unit, f, t) }
        + v_startup_MIP(unit, starttype, s, f, t)
            ${ uft_onlineMIP(unit, f, t) }
      ) // END sum(starttype)
        * p_gnu(grid, node, unit, 'unitSizeTot')
        * (
            + sum(suft(effGroup, unit, f, t), // Uses the minimum 'lb' for the current efficiency approximation