Making it possible to choose in the modelsInit.gms file whether to include the...

Making it possible to choose in the modelsInit.gms file whether to include the two additional incremental heat rate equations that use binary variables:
- mSettings(m, 'incHRAdditionalConstraints') = 0: use the constraints but only for units with non-convex fuel use
- mSettings(m, 'incHRAdditionalConstraints') = 1: use the constraints for all units represented using incremental heat rates
Enabling incremental heat rates also in the invest model (and potentially in the building model). Some clean-up to the incremental heat rate equations.

defModels/building.gms

@@ -45,6 +45,11 @@ Model building /
 *    q_conversionSOS2InputIntermediate
 *    q_conversionSOS2Constraint
 *    q_conversionSOS2IntermediateOutput
+*    q_conversionIncHR
+*    q_conversionIncHRMaxGen
+*    q_conversionIncHRBounds
+*    q_conversionIncHR_help1
+*    q_conversionIncHR_help2
 *    q_fuelUseLimit
 
     // Energy Transfer

defModels/buildingInit_temp.gms

@@ -117,6 +117,13 @@ if (mType('building'),
     // Define the length of the initialization period. Results outputting starts after the period. Uses ord(t) > t_start + t_initializationPeriod in the code.
     mSettings('building', 't_initializationPeriod') = 0;  // r_state and r_online are stored also for the last step in the initialization period, i.e. ord(t) = t_start + t_initializationPeriod
 
+* --- Define the use of additional constraints for units with incremental heat rates
+
+    // How to use q_conversionIncHR_help1 and q_conversionIncHR_help2
+    mSettings('building', 'incHRAdditionalConstraints') = 0;
+    // 0 = use the constraints but only for units with non-convex fuel use
+    // 1 = use the constraints for all units represented using incremental heat rates
+
 * --- Control the solver ------------------------------------------------------
 
     // Control the use of advanced basis

defModels/invest.gms

@@ -47,6 +47,11 @@ Model invest /
     q_conversionSOS2InputIntermediate
     q_conversionSOS2Constraint
     q_conversionSOS2IntermediateOutput
+    q_conversionIncHR
+    q_conversionIncHRMaxGen
+    q_conversionIncHRBounds
+    q_conversionIncHR_help1
+    q_conversionIncHR_help2
     q_fuelUseLimit
 
     // Energy Transfer

defModels/investInit_temp.gms

@@ -123,6 +123,13 @@ if (mType('invest'),
     // Define the length of the initialization period. Results outputting starts after the period. Uses ord(t) > t_start + t_initializationPeriod in the code.
     mSettings('invest', 't_initializationPeriod') = 0;  // r_state and r_online are stored also for the last step in the initialization period, i.e. ord(t) = t_start + t_initializationPeriod
 
+* --- Define the use of additional constraints for units with incremental heat rates
+
+    // How to use q_conversionIncHR_help1 and q_conversionIncHR_help2
+    mSettings('invest', 'incHRAdditionalConstraints') = 0;
+    // 0 = use the constraints but only for units with non-convex fuel use
+    // 1 = use the constraints for all units represented using incremental heat rates
+
 * --- Control the solver ------------------------------------------------------
 
     // Control the use of advanced basis

defModels/schedule.gms

@@ -50,8 +50,8 @@ Model schedule /
     q_conversionIncHR
     q_conversionIncHRMaxGen
     q_conversionIncHRBounds
-*    q_conversionIncHR_help1
-*    q_conversionIncHR_help2
+    q_conversionIncHR_help1
+    q_conversionIncHR_help2
     q_fuelUseLimit
 
     // Energy Transfer

defModels/scheduleInit_temp.gms

@@ -172,6 +172,13 @@ if (mType('schedule'),
     // Define the length of the initialization period. Results outputting starts after the period. Uses ord(t) > t_start + t_initializationPeriod in the code.
     mSettings('schedule', 't_initializationPeriod') = 0;  // r_state and r_online are stored also for the last step in the initialization period, i.e. ord(t) = t_start + t_initializationPeriod
 
+* --- Define the use of additional constraints for units with incremental heat rates
+
+    // How to use q_conversionIncHR_help1 and q_conversionIncHR_help2
+    mSettings('schedule', 'incHRAdditionalConstraints') = 0;
+    // 0 = use the constraints but only for units with non-convex fuel use
+    // 1 = use the constraints for all units represented using incremental heat rates
+
 * --- Control the solver ------------------------------------------------------
 
     // Control the use of advanced basis

defOutput/debugSymbols.inc

@@ -91,6 +91,7 @@ v_investTransfer_LP
 v_investTransfer_MIP
 v_invest_LP
 v_invest_MIP
+v_help_inc
 
 * equations
     q_obj
@@ -118,6 +119,11 @@ v_invest_MIP
     q_conversionSOS2InputIntermediate
     q_conversionSOS2Constraint
     q_conversionSOS2IntermediateOutput
+    q_conversionIncHR
+    q_conversionIncHRMaxGen
+    q_conversionIncHRBounds
+    q_conversionIncHR_help1
+    q_conversionIncHR_help2
     q_fuelUseLimit
 
     // Energy Transfer

inc/1a_definitions.gms

@@ -58,6 +58,8 @@ Sets
         dataLength, // The maximum number of time steps in any input data time series (recommended for correctly circulating data)
         red_num_leaves "Desired number of preserved scenarios or leaves (SCENRED)"
         red_percentage "Desired relative distance (accuracy) of scenario reduction (SCENRED)"
+        incHRAdditionalConstraints // Method to include the two additional constraints for incremental heat rates;
+                                   // 0 = include for units with non-convex fuel use, 1 = include for all units
         /
 
     // Solve info
@@ -297,7 +299,7 @@ param_unit "Set of possible data parameters for units" /
     eff00 * eff12 "Efficiency of the unit to convert input to output/intermediate product"
     opFirstCross  "The operating point where the real efficiency curve and approximated efficiency curve cross"
     op00 * op12   "Right border of the efficiency point"
-    hr00 * hr12   "incremental heat rates"
+    hr00 * hr12   "Incremental heat rates (GJ/MWh)"
     hrop00 * hrop12   "Right border of the incremental heat rates"
     section       "Possibility to define a no load fuel use for units with zero minimum output"
     hrsection     "no load fuel use to be defined when using incremental heat rates"

inc/1b_sets.gms

@@ -49,6 +49,7 @@ Sets
     unit_investLP(unit) "Units with continuous investments allowed"
     unit_investMIP(unit) "Units with integer investments allowed"
     unit_timeseries(unit) "Units with time series enabled"
+    unit_incHRAdditionalConstraints(unit) "Units that use the two additional incremental heat rate constraints" 
 
 * --- Nodes -------------------------------------------------------------------
     node_spill(node) "Nodes that can spill; used to remove v_spill variables where not relevant"

inc/2b_eqDeclarations.gms

@@ -23,10 +23,12 @@ $offtext
 
 Scalars
     PENALTY "Default equation violation penalty"
+    BIG_M "A large number used together with with binary variables in some equations"
 ;
 
 $If set penalty PENALTY=%penalty%;
 $If not set penalty PENALTY=1e9;
+BIG_M = 1e5;
 
 Parameters
     PENALTY_BALANCE(grid, node) "Penalty on violating energy balance eq. (EUR/MWh)"
@@ -78,9 +80,6 @@ equations
     q_conversionIncHR(s, effSelector, unit, f, t)  "Conversion of inputs to outputs for incremental heat rates"
     q_conversionIncHRMaxGen(grid, node,s, effSelector, unit, f, t)  "Max Generating level"
     q_conversionIncHRBounds(grid, node, s, hr, effSelector, unit, f, t) "Heat rate bounds"
-    q_conversionIncHRBounds1(grid, node,s, effSelector, unit, f, t) "Heat rate bounds for hr=1"
-    q_conversionIncHRBounds2(grid, node, s, effSelector, unit, f, t) "Heat rate bounds for hr=2"
-    q_conversionIncHRBounds3(grid, node, s, effSelector, unit, f, t) "Heat rate bounds for hr=3"
     q_conversionIncHR_help1(grid, node, s, hr, effSelector, unit, f, t) "Helper equation 1 to ensure that the first heat rate segments are used first"
     q_conversionIncHR_help2(grid, node, s, hr, effSelector, unit, f, t) "Helper equation 2 to ensure that the first heat rate segments are used first"
     q_fuelUseLimit(s, fuel, unit, f, t) "Fuel use cannot exceed limits"

inc/2d_constraints.gms

@@ -950,7 +950,6 @@ q_conversionDirectInputOutput(s, suft(effDirect(effGroup), unit, f, t))$sft(s, f
 ;
 * --- Incremental Heat Rate Conversion ------------------------------------------
 
-
 q_conversionIncHR(s, suft(effIncHR(effGroup), unit, f, t))$sft(s, f, t) ..
 
     // Sum over endogenous energy inputs
@@ -967,21 +966,20 @@ q_conversionIncHR(s, suft(effIncHR(effGroup), unit, f, t))$sft(s, f, t) ..
 
     // Sum over energy outputs
     + sum(gnu_output(grid, node, unit),
-        + sum (hr,
-         v_gen_inc(grid, node, unit, hr, s, f, t) //   output of each heat rate segment
-
-
-            * [ // heat rate
-                p_unit(unit, hr)/3.6     // p/ts?
+        + sum(hr,
+            + v_gen_inc(grid, node, unit, hr, s, f, t) // output of each heat rate segment
+            * [
+                + p_unit(unit, hr) // heat rate
+                / 3.6 // unit conversion from [GJ/MWh] into [MWh/MWh]
                 ] // END * v_gen_inc
-              )) // END sum(gnu_output)
+            ) // END sum(hr)
+        ) // END sum(gnu_output)
 
     // Consumption of keeping units online (no-load fuel use)
     + sum(gnu_output(grid, node, unit),
         + p_gnu(grid, node, unit, 'unitSizeGen')
         ) // END sum(gnu_output)
         * [ // Unit online state
-            + 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)}
 
             // Run-up and shutdown phase efficiency correction
@@ -1007,21 +1005,23 @@ q_conversionIncHR(s, suft(effIncHR(effGroup), unit, f, t))$sft(s, f, t) ..
 ;
 
 * --- Incremental Heat Rate Conversion ------------------------------------------
+
 q_conversionIncHRMaxGen(gn(grid, node), s, suft(effIncHR(effGroup), unit, f, t))${  sft(s, f, t)
                                                                                     and gnu_output(grid, node, unit)
                                                                                     } ..
 
-    v_gen(grid, node, unit, s, f, t)
+    + v_gen(grid, node, unit, s, f, t)
 
     =E=
 
     // Sum over heat rate segments
-    sum(hr$(p_unit(unit, hr)),
-      v_gen_inc(grid, node, unit, hr, s, f, t)
-      )// END sum (hr)
+    + sum(hr$(p_unit(unit, hr)),
+        + v_gen_inc(grid, node, unit, hr, s, f, t)
+        )// END sum (hr)
 ;
 
 * --- Incremental Heat Rate Conversion ------------------------------------------
+
 q_conversionIncHRBounds(gn(grid, node), s, hr, suft(effIncHR(effGroup), unit, f, t))${  sft(s, f, t)
                                                                                         and gnu_output(grid, node, unit)
                                                                                         and p_unit(unit, hr)
@@ -1039,26 +1039,14 @@ q_conversionIncHRBounds(gn(grid, node), s, hr, suft(effIncHR(effGroup), unit, f,
         *  v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
 ;
 
-* --- Incremental Heat Rate Conversion ------------------------------------------
-q_conversionIncHRBounds1(grid, node,s,  suft(effIncHR(effGroup), unit, f, t))$(sft(s, f, t)) ..
-
-    v_gen_inc(grid, node, unit, 'hr00', s, f, t)
-
-    =L=
-    (
-      p_effUnit(effGroup, unit, effGroup, 'lb')
-      )
-       *  p_gnu(grid, node, unit, 'maxGen')*v_online_MIP(unit, s, f+df_central(f,t), t)${uft_onlineMIP(unit, f, t)}
-;
-
-
 * --- Incremental Heat Rate Conversion (First Segments First) -----------------
 
-q_conversionIncHR_help1(gn(grid, node), s, hr, suft(effIncHR(effGroup), unit, f, t))${  sft(s, f, t)
-                                                                                        and gnu_output(grid, node, unit)
-                                                                                        and p_unit(unit, hr)
-                                                                                        and p_unit(unit, hr+1)
-                                                                                        } ..
+q_conversionIncHR_help1(gn(grid, node), s, hr, suft(effIncHR(effGroup), unit_incHRAdditionalConstraints(unit), f, t))
+    ${  sft(s, f, t)
+        and gnu_output(grid, node, unit)
+        and p_unit(unit, hr)
+        and p_unit(unit, hr+1)
+        } ..
 
     + v_gen_inc(grid, node, unit, hr, s, f, t)
     - (
@@ -1070,21 +1058,22 @@ q_conversionIncHR_help1(gn(grid, node), s, hr, suft(effIncHR(effGroup), unit, f,
 
     =G=
 
-    - 100000  // a very large number, add to scalars in 2b_eqDeclarations.gms?
+    - BIG_M
         * (1 - v_help_inc(grid, node, unit, hr, s, f, t))
 ;
 
-q_conversionIncHR_help2(gn(grid, node), s, hr, suft(effIncHR(effGroup), unit, f, t))${  sft(s, f, t)
-                                                                                        and gnu_output(grid, node, unit)
-                                                                                        and p_unit(unit, hr)
-                                                                                        and p_unit(unit, hr-1)
-                                                                                        } ..
+q_conversionIncHR_help2(gn(grid, node), s, hr, suft(effIncHR(effGroup), unit_incHRAdditionalConstraints(unit), f, t))
+    ${  sft(s, f, t)
+        and gnu_output(grid, node, unit)
+        and p_unit(unit, hr)
+        and p_unit(unit, hr-1)
+        } ..
 
     + v_gen_inc(grid, node, unit, hr, s, f, t)
 
     =L=
 
-    + 100000  // a very large number, add to scalars in 2b_eqDeclarations.gms?
+    + BIG_M
         * v_help_inc(grid, node, unit, hr-1, s, f, t)
 ;
 

inc/3a_periodicInit.gms

@@ -411,6 +411,19 @@ loop(effGroupSelectorUnit(effSelector, unit, effSelector_),
         p_effUnit(effSelector, unit, effSelector, 'op') = smax(hrop, p_unit(unit, hrop)); // Maximum operating point
         p_effUnit(effSelector, unit, effSelector, 'slope') = 1 / smax(eff${p_unit(unit, eff)}, p_unit(unit, eff)); // Uses maximum found (nonzero) efficiency.
         p_effUnit(effSelector, unit, effSelector, 'section') = p_unit(unit, 'hrsection'); // pre-defined
+
+        // Whether to use q_conversionIncHR_help1 and q_conversionIncHR_help2 or not
+        loop(m,
+            loop(hr${p_unit(unit, hr)},
+                if (mSettings(m, 'incHRAdditionalConstraints') = 0,
+                    if (p_unit(unit, hr) < p_unit(unit, hr-1),
+                        unit_incHRAdditionalConstraints(unit) = yes;
+                    ); // END if(hr)
+                else
+                    unit_incHRAdditionalConstraints(unit) = yes;
+                ); // END if(incHRAdditionalConstraints)
+            ); // END loop(hr)
+        ); // END loop(m)
     ); // END if(effIncHR)
 ); // END loop(effGroupSelectorUnit)
 

inc/3b_periodicLoop.gms

@@ -33,6 +33,8 @@ $ifthene.debug not %debug%>0
     Option clear = v_online_MIP;
     Option clear = v_invest_MIP;
     Option clear = v_investTransfer_MIP;
+    // Binary Variables
+    Option clear = v_help_inc;
     // SOS2 Variables
     Option clear = v_sos2;
     // Positive Variables
@@ -49,6 +51,7 @@ $ifthene.debug not %debug%>0
     Option clear = v_investTransfer_LP;
     Option clear = v_online_LP;
     Option clear = v_invest_LP;
+    Option clear = v_gen_inc;
     // Feasibility control
     Option clear = v_stateSlack;
     Option clear = vq_gen;
@@ -84,6 +87,11 @@ $ifthene.debug not %debug%>0
     Option clear = q_conversionSOS2InputIntermediate;
     Option clear = q_conversionSOS2Constraint;
     Option clear = q_conversionSOS2IntermediateOutput;
+    Option clear = q_conversionIncHR;
+    Option clear = q_conversionIncHRMaxGen;
+    Option clear = q_conversionIncHRBounds;
+    Option clear = q_conversionIncHR_help1;
+    Option clear = q_conversionIncHR_help2;
     Option clear = q_fuelUseLimit;
 
     // Energy Transfer