Merge pull request #772 from spine-tools/min_res_gen_constraint_as_cuts

Turn min-res-gen constraint into cuts and report constraint value

src/SpineOpt.jl

@@ -187,7 +187,7 @@ include("constraints/constraint_min_node_voltage_angle.jl")
 include("constraints/constraint_node_voltage_angle.jl")
 include("constraints/constraint_connection_unitary_gas_flow.jl")
 include("constraints/constraint_mp_any_invested_cuts.jl")
-include("constraints/constraint_mp_min_res_gen_to_demand_ratio.jl")
+include("constraints/constraint_mp_min_res_gen_to_demand_ratio_cuts.jl")
 include("constraints/constraint_investment_group_equal_investments.jl")
 include("constraints/constraint_investment_group_entities_invested_available.jl")
 include("constraints/constraint_investment_group_capacity_invested_available.jl")

src/constraints/constraint_mp_any_invested_cuts.jl

@@ -33,7 +33,7 @@ function add_constraint_mp_any_invested_cuts!(m::Model)
     ) = m.ext[:spineopt].variables
     merge!(
         get!(m.ext[:spineopt].constraints, :mp_any_invested_cut, Dict()),
-            Dict(
+        Dict(
             (benders_iteration=bi, t=t) => @constraint(
                 m,
                 + sp_objective_upperbound[t]

src/constraints/constraint_mp_min_res_gen_to_demand_ratio_cuts.jl

@@ -0,0 +1,95 @@
+#############################################################################
+# Copyright (C) 2017 - 2023  Spine Project
+#
+# This file is part of SpineOpt.
+#
+# SpineOpt 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.
+#
+# SpineOpt 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+#############################################################################
+
+"""
+    add_constraint_mp_min_res_gen_to_demand_ratio!(m::Model)
+
+sum (
+    + res generation from subproblem
+    + (units_invested_available - units_invested_available from last iteration)
+    * unit_availability_factor * unit_capacity * unit_conv_cap_to_flow
+) >= mp_min_res_gen_to_demand_ratio * total demand
+"""
+function add_constraint_mp_min_res_gen_to_demand_ratio_cuts!(m::Model)
+    @fetch units_invested_available, mp_min_res_gen_to_demand_ratio_slack = m.ext[:spineopt].variables
+    merge!(
+        get!(m.ext[:spineopt].constraints, :mp_min_res_gen_to_demand_ratio_cuts, Dict()),
+        Dict(
+            (benders_iteration=bi, commodity=comm) => @constraint(
+                m,
+                + sum(
+                    window_sum_duration(m, sp_unit_flow(unit=u, node=n, direction=d, stochastic_scenario=s), window)
+                    for window in m.ext[:spineopt].temporal_structure[:sp_windows]
+                    for (u, s) in unit_stochastic_indices(m; unit=unit(is_renewable=true))
+                    for (u, n, d) in unit__to_node(unit=u, node=node__commodity(commodity=comm), _compact=false);
+                    init=0
+                )
+                + sum(
+                    sum(
+                        + units_invested_available[u, s, t]
+                        - internal_fix_units_invested_available(unit=u, stochastic_scenario=s, t=t, _default=0)
+                        for (u, s, t) in units_invested_available_indices(
+                            m; unit=u, stochastic_scenario=s, t=to_time_slice(m; t=window)
+                        );
+                        init=0
+                    )
+                    * window_sum_duration(
+                        m,
+                        + unit_availability_factor(unit=u, stochastic_scenario=s)
+                        * unit_capacity(unit=u, node=n, direction=d, stochastic_scenario=s)
+                        * unit_conv_cap_to_flow(unit=u, node=n, direction=d, stochastic_scenario=s),
+                        window
+                    )
+                    for window in m.ext[:spineopt].temporal_structure[:sp_windows]
+                    for (u, s) in unit_stochastic_indices(m; unit=unit(is_renewable=true)) 
+                    for (u, n, d) in unit__to_node(unit=u, node=node__commodity(commodity=comm), _compact=false);
+                    init=0,
+                )
+                + get(mp_min_res_gen_to_demand_ratio_slack, (comm,), 0)
+                >=
+                + mp_min_res_gen_to_demand_ratio(commodity=comm)
+                * (
+                    sum(
+                        window_sum_duration(m, demand(node=n, stochastic_scenario=s), window)
+                        for window in m.ext[:spineopt].temporal_structure[:sp_windows]
+                        for (n, s) in node_stochastic_indices(
+                            m; node=intersect(indices(demand), node__commodity(commodity=comm))
+                        );
+                        init=0
+                    )
+                    + sum(
+                        window_sum_duration(
+                            m,
+                            fractional_demand(node=n, stochastic_scenario=s) * demand(node=ng, stochastic_scenario=s),
+                            window
+                        )
+                        for window in m.ext[:spineopt].temporal_structure[:sp_windows]
+                        for (n, s) in node_stochastic_indices(
+                            m; node=intersect(indices(fractional_demand), node__commodity(commodity=comm))
+                        )
+                        for ng in groups(n);
+                        init=0
+                    )
+                )
+            )
+            for bi in last(benders_iteration())
+            for comm in indices(mp_min_res_gen_to_demand_ratio)
+        )
+    )
+end

src/run_spineopt_benders.jl

@@ -101,9 +101,6 @@ function rerun_spineopt_benders!(
             break
         end
         @timelog log_level 2 "Add MP cuts..." _add_mp_cuts!(m_mp; log_level=log_level)
-        @timelog log_level 2 "Adding MP renewing constraints...\n" _add_mp_renewing_constraints!(
-            m_mp; log_level=log_level
-        )
         _unfix_history!(m)
         j += 1
         global current_bi = add_benders_iteration(j)
@@ -172,17 +169,6 @@ function _add_mp_constraints!(m; log_level=3)
     _update_constraint_names!(m)
 end
 
-function _add_mp_renewing_constraints!(m; log_level=3)
-    for add_constraint! in (
-            add_constraint_mp_min_res_gen_to_demand_ratio!,
-        )
-        name = name_from_fn(add_constraint!)
-        @timelog log_level 3 "- [$name]" add_constraint!(m)
-    end
-    _update_constraint_names!(m)
-end
-
-
 function _add_constraint_sp_objective_upperbound!(m::Model)
     @fetch sp_objective_upperbound = m.ext[:spineopt].variables
     m.ext[:spineopt].constraints[:mp_objective] = Dict(
@@ -218,12 +204,14 @@ end
 Add benders cuts to master problem.
 """
 function _add_mp_cuts!(m; log_level=3)
-    @timelog log_level 3 " - [constraint_mp_any_invested_cuts]" add_constraint_mp_any_invested_cuts!(m)
-    # Name constraints
-    cons = m.ext[:spineopt].constraints[:mp_any_invested_cut]
-    for (inds, con) in cons
-        _set_name(con, string(:mp_any_invested_cut, inds))
+    for add_constraint! in (
+            add_constraint_mp_any_invested_cuts!,
+            add_constraint_mp_min_res_gen_to_demand_ratio_cuts!,
+        )
+        name = name_from_fn(add_constraint!)
+        @timelog log_level 3 "- [$name]" add_constraint!(m)
     end
+    _update_constraint_names!(m)
 end
 
 function _unfix_history!(m::Model)

src/run_spineopt_mga.jl

@@ -94,20 +94,10 @@ function rerun_spineopt_mga!(
             && isempty(indices(storages_invested_big_m_mga))
             && mga_iterations < max_mga_iteration
         )
-            for cons in
-                [:mga_objective_ub,
-                :mga_diff_ub1,]
-                for k in keys(m.ext[:spineopt].constraints[cons])
-                    try m.ext[:spineopt].constraints[cons][k]
-                        delete(m, m.ext[:spineopt].constraints[cons][k])
-                    catch
-                    end
-                end
-            end
-            for vars in  [:mga_aux_diff,]
-                for k in keys(m.ext[:spineopt].variables[vars])
-                    try m.ext[:spineopt].constraints[cons][k]
-                        delete(m, m.ext[:spineopt].variables[vars][k])
+            for name in (:mga_objective_ub, :mga_diff_ub1)
+                for con in values(m.ext[:spineopt].constraints[name])
+                    try
+                        delete(m, con)
                     catch
                     end
                 end

src/run_spineopt_standard.jl

@@ -270,7 +270,7 @@ function run_spineopt_kernel!(
     k = _resume_run!(m, resume_file_path, log_level, update_names)
     k === nothing && return m
     calculate_duals = any(
-        startswith(lowercase(name), r"bound_|constraint_") for name in String.(keys(m.ext[:spineopt].outputs))
+        startswith(name, r"bound_|constraint_") for name in lowercase.(string.(keys(m.ext[:spineopt].outputs)))
     )
     while optimize
         @log log_level 1 "\nWindow $k: $(current_window(m))"
@@ -381,6 +381,7 @@ Save a model results: first postprocess results, then save variables and objecti
 """
 function _save_model_results!(m; iterations=nothing)
     _save_variable_values!(m)
+    _save_constraint_values!(m)
     _save_objective_values!(m)
 end
 
@@ -393,6 +394,17 @@ function _save_variable_values!(m::Model)
     end
 end
 
+"""
+Save the value of all constraints if the user wants to report it.
+"""
+function _save_constraint_values!(m::Model)
+    for (name, con) in m.ext[:spineopt].constraints
+        name = Symbol(:value_constraint_, name)
+        name in keys(m.ext[:spineopt].outputs) || continue
+        m.ext[:spineopt].values[name] = Dict(ind => JuMP.value(c) for (ind, c) in con)
+    end
+end
+
 """
 The value of a JuMP variable, rounded if necessary.
 """
@@ -916,10 +928,21 @@ function update_model!(m; log_level=3, update_names=false)
     @timelog log_level 2 "Applying non-anticipativity constraints..." apply_non_anticipativity_constraints!(m)
 end
 
-function _update_constraint_names!(m)
-    for (key, cons) in m.ext[:spineopt].constraints                        
-        for (ind, con) in cons        
-            constraint_name = _sanitize_constraint_name(string(key, ind))                            
+function _update_variable_names!(m, names=keys(m.ext[:spineopt].variables))
+    for name in names   
+        var = m.ext[:spineopt].variables[name]
+        # NOTE: only update names for the representative variables
+        # This is achieved by using the indices function from the variable definition
+        for ind in m.ext[:spineopt].variables_definition[name][:indices](m)
+            _set_name(var[ind], _base_name(name, ind))
+        end
+    end
+end
+
+function _update_constraint_names!(m, names=keys(m.ext[:spineopt].constraints))
+    for name in names   
+        for (ind, con) in m.ext[:spineopt].constraints[name]        
+            constraint_name = _sanitize_constraint_name(string(name, ind))                            
             _set_name(con, constraint_name)
         end
     end
@@ -931,16 +954,6 @@ function _sanitize_constraint_name(constraint_name)
     replace(constraint_name, pattern => "_")
 end
 
-function _update_variable_names!(m)
-    for (name, var) in m.ext[:spineopt].variables
-        # NOTE: only update names for the representative variables
-        # This is achieved by using the indices function from the variable definition
-        for ind in m.ext[:spineopt].variables_definition[name][:indices](m)
-            _set_name(var[ind], _base_name(name, ind))
-        end
-    end
-end
-
 _set_name(x::Union{VariableRef,ConstraintRef}, name) = set_name(x, name)
 _set_name(::Union{Call,Nothing}, name) = nothing