Fix PRICE_EMISSION calculation bug

message_ix/model/MESSAGE/data_load.gms

@@ -97,6 +97,8 @@ Set rating_unrated(rating) ;
 rating_unrated(rating) = yes ;
 rating_unrated('unrated') = no ;
 
+Set newtec(tec) / wind_ppl / ;
+
 *----------------------------------------------------------------------------------------------------------------------*
 * assignment and computation of MESSAGE-specific auxiliary parameters                                                  *
 *----------------------------------------------------------------------------------------------------------------------*

message_ix/model/MESSAGE/model_learningeos.gms

@@ -0,0 +1,47 @@
+SETS
+  size           'size'  / small, medium, large / ;
+
+ALIAS (size,size2);
+PARAMETERS
+  cap_new2(node,newtec,year_all2)        'annual newly installed capacity'
+  bin_cap_new(node,newtec,year_all2)     'binary of newly installed capacity'
+  rho(newtec)                            'economy of scale parameter'            / wind_ppl      1.0       / #0.8
+  b(newtec)                              'technology cost learning parameter'    / wind_ppl      0.9     / #0.9
+  u(size)                                'unit size'
+         / small      5
+           medium     10
+           large      50     /
+  inv_cost_ref(node,newtec)              'initial capex'
+  nbr_unit_ref(newtec)                   'initial number of unit'                / wind_ppl      100     /
+  u_ref(newtec)                          'reference size'                        / wind_ppl      5       / ;
+inv_cost_ref(node,newtec) = 1500;
+
+SCALAR hist_length                       the length of historical periods;
+hist_length = card(year_all2) - card(model_horizon);
+
+VARIABLES
+  NBR_UNIT(node,newtec,size,year_all2)   number of units for each size every year
+  CAPEX_TEC(node,newtec,year_all2)       capital cost in dollar per kW
+  OBJECT                                 objective function ;
+
+POSITIVE VARIABLES
+  NBR_UNIT ;
+
+EQUATIONS
+  OBJECTIVE_INNER        total investment cost
+  CAP_NEW_BALANCE        installed capacity balance
+  CAPEX_ESTIMATE         estimating average capex
+  NO_BUILT_YEAR          annual investment cost
+;
+
+
+OBJECTIVE_INNER..                        OBJECT =e= sum((node,newtec,year_all2), CAPEX_TEC(node,newtec,year_all2)*cap_new2(node,newtec,year_all2)) ;
+CAP_NEW_BALANCE(node,newtec,year_all2).. sum(size, NBR_UNIT(node,newtec,size,year_all2)*u(size)) =e= cap_new2(node,newtec,year_all2) ;
+CAPEX_ESTIMATE(node,newtec,year_all2)..  CAPEX_TEC(node,newtec,year_all2)*cap_new2(node,newtec,year_all2) =g= sum(size,inv_cost_ref(node,newtec)
+                                              * NBR_UNIT(node,newtec,size,year_all2)*u(size)
+                                              * [(((sum((size2,year_all3)$(ord(year_all3) le ord(year_all2) and ord(year_all3) gt hist_length), NBR_UNIT(node,newtec,size2,year_all3))+nbr_unit_ref(newtec))/nbr_unit_ref(newtec))**(-b(newtec)))]
+                                              * [((u(size)/u_ref(newtec))**rho(newtec))]) ;
+NO_BUILT_YEAR(node,newtec,year_all2)..   CAPEX_TEC(node,newtec,year_all2) =e= bin_cap_new(node,newtec,year_all2)*CAPEX_TEC(node,newtec,year_all2)
+                                                                              + (1-bin_cap_new(node,newtec,year_all2))*CAPEX_TEC(node,newtec,year_all2-1) ;
+
+model learningeos / all /;

message_ix/model/MESSAGE/model_setup.gms

@@ -21,6 +21,11 @@ $IF NOT SET out      $SETGLOBAL out "output/MsgOutput.gdx"
 * rolling horizon (period-by-period, recursive-dynamic with limited foresight - 'number of years of foresight'
 $IF NOT SET foresight   $SETGLOBAL foresight "0"
 
+** define learning mode (active / inactive) **
+* deactivate - 0 (assumed as default if not specified)
+* activate - 1
+$IF NOT SET learningmode   $SETGLOBAL learningmode "0"
+
 ** specify optional additional calibration output **
 $IF NOT SET calibration $SETGLOBAL calibration ""
 * mark with * to include detailed calibration information in outputs and get an extended GAMS listing (.lst) file
@@ -84,3 +89,4 @@ $INCLUDE MESSAGE/scaling_investment_costs.gms
 *----------------------------------------------------------------------------------------------------------------------*
 
 $INCLUDE MESSAGE/model_core.gms
+$INCLUDE MESSAGE/model_learningeos.gms

message_ix/model/MESSAGE/model_solve.gms

@@ -5,6 +5,11 @@
 * This part of the code includes the perfect-foresight, myopic and rolling-horizon model solve statements
 * including the required accounting of investment costs beyond the model horizon.
 ***
+*set year_rd(year_all) /700, 710/;
+*set year_rd(year_all) all year in recursive dynamic iterations ;
+*year_rd(year_all) = ord(year_all) >1;
+Parameter iter 'iteration number';
+iter = 1;
 
 if (%foresight% = 0,
 ***
@@ -24,6 +29,7 @@ if (%foresight% = 0,
 
 * write a status update to the log file, solve the model
     put_utility 'log' /'+++ Solve the perfect-foresight version of MESSAGEix +++ ' ;
+    option threads = 4 ;
     Solve MESSAGE_LP using LP minimizing OBJ ;
 
 * write model status summary
@@ -55,7 +61,7 @@ EMISSION_CONSTRAINT.m(node,type_emission,type_tec,type_year)$(
     PRICE_EMISSION.l(node,type_emission,type_tec,year)$( SUM(type_year$( cat_year(type_year,year) ), 1 ) ) =
         SMAX(type_year$( cat_year(type_year,year) ),
                - EMISSION_CONSTRAINT.m(node,type_emission,type_tec,type_year) )
-            / df_year(year) ;
+             * (sum(year_all$(ord(year_all) = 1), df_period(year_all))/df_period(year)) * (duration_period(year)/(sum(year_all$(ord(year_all) = 1), duration_period(year_all)))) ;
     PRICE_EMISSION.l(node,type_emission,type_tec,year)$(
         PRICE_EMISSION.l(node,type_emission,type_tec,year) = - inf ) = 0 ;
 
@@ -87,47 +93,56 @@ else
 * variables in additional reporting parameters - the last model solve automatically includes the results over the
 * entire model horizon and can be imported via the ixmp interface.
 ***
-
     year(year_all) = no ;
 
     LOOP(year_all$( model_horizon(year_all) ),
 
 * include all past periods and future periods including the period where the %foresight% is reached
-        year(year_all) = yes ;
-
-* reset the investment cost scaling parameter
-        year(year_all2)$( ORD(year_all2) > ORD(year_all)
-            AND duration_period_sum(year_all,year_all2) < %foresight% ) = yes ;
-
-* write a status update and time elapsed to the log file, solve the model
-        put_utility 'log' /'+++ Solve the recursive-dynamic version of MESSAGEix - iteration ' year_all.tl:0 '  +++ ' ;
-        $$INCLUDE includes/aux_computation_time.gms
-        Solve MESSAGE_LP using LP minimizing OBJ ;
-
-* write model status summary
-        status(year_all,'modelstat') = MESSAGE_LP.modelstat ;
-        status(year_all,'solvestat') = MESSAGE_LP.solvestat ;
-        status(year_all,'resUsd')    = MESSAGE_LP.resUsd ;
-        status(year_all,'objEst')    = MESSAGE_LP.objEst ;
-        status(year_all,'objVal')    = MESSAGE_LP.objVal ;
-
-* write an error message AND ABORT THE SOLVE LOOP if model did not solve to optimality
-        IF( NOT ( MESSAGE_LP.modelstat = 1 OR MESSAGE_LP.modelstat = 8 ),
-            put_utility 'log' /'+++ MESSAGEix did not solve to optimality - run is aborted, no output produced! +++ ' ;
-            ABORT "MESSAGEix did not solve to optimality!"
-        ) ;
+             year(year_all2)$( ORD(year_all2) < (ORD(year_all) + %foresight%) ) = yes ;
+             year4(year_all2)$((ord(year_all2) < ord(year_all))) = yes ;
+
+             option threads = 4 ;
+             Solve MESSAGE_LP using LP minimizing OBJ ;
+*            write model status summary
+             status('perfect_foresight','modelstat') = MESSAGE_LP.modelstat ;
+             status('perfect_foresight','solvestat') = MESSAGE_LP.solvestat ;
+             status('perfect_foresight','resUsd')    = MESSAGE_LP.resUsd ;
+             status('perfect_foresight','objEst')    = MESSAGE_LP.objEst ;
+             status('perfect_foresight','objVal')    = MESSAGE_LP.objVal ;
+
+*            write an error message if model did not solve to optimality
+             IF( NOT ( MESSAGE_LP.modelstat = 1 OR MESSAGE_LP.modelstat = 8 ),
+                 put_utility 'log' /'+++ MESSAGEix did not solve to optimality - run is aborted, no output produced! +++ ' ;
+                 ABORT "MESSAGEix did not solve to optimality!"
+             ) ;
+
+             IF(%learningmode% = 1,
+*            passing CAP_NEW values to update cap_new2 data for unit and size optimization
+                 cap_new2(node,newtec,year_all2) = CAP_NEW.l(node,newtec,year_all2);
+*            this is to make bin param equal to 1 when technology is built, and 0 if otherwise
+                 bin_cap_new(node,newtec,year_all2) = CAP_NEW.l(node,newtec,year_all2);
+                 bin_cap_new(node,newtec,year_all2)$(bin_cap_new(node,newtec,year_all2) > 0) = 1 ;
+                 solve learningeos using nlp minimizing OBJECT;
+*            passing CapexTec values to update inv_cost data for MESSAGE optimization
+                 inv_cost(node,newtec,year_all2) = CAPEX_TEC.l(node,newtec,year_all2);
+
+                 display bin_cap_new, NBR_UNIT.l, CAPEX_TEC.l, inv_cost, cap_new2, CAP_NEW.l;
+                 );
 
 * fix all variables of the current iteration period 'year_all' to the optimal levels
-        EXT.fx(node,commodity,grade,year_all) =  EXT.l(node,commodity,grade,year_all) ;
-        CAP_NEW.fx(node,tec,year_all) = CAP_NEW.l(node,tec,year_all) ;
-        CAP.fx(node,tec,year_all2,year_all)$( map_period(year_all2,year_all) ) = CAP.l(node,tec,year_all,year_all2) ;
-        ACT.fx(node,tec,year_all2,year_all,mode,time)$( map_period(year_all2,year_all) )
-            = ACT.l(node,tec,year_all2,year_all,mode,time) ;
-        CAP_NEW_UP.fx(node,tec,year_all) = CAP_NEW_UP.l(node,tec,year_all) ;
-        CAP_NEW_LO.fx(node,tec,year_all) = CAP_NEW_LO.l(node,tec,year_all) ;
-        ACT_UP.fx(node,tec,year_all,time) = ACT_UP.l(node,tec,year_all,time) ;
-        ACT_LO.fx(node,tec,year_all,time) = ACT_LO.l(node,tec,year_all,time) ;
-
+        EXT.fx(node,commodity,grade,year4) =  EXT.l(node,commodity,grade,year4) ;
+        CAP_NEW.fx(node,tec,year4) = CAP_NEW.l(node,tec,year4) ;
+*        CAP.fx(node,tec,year4,year4) = CAP.l(node,tec,year4,year4) ;
+        CAP.up(node,tec,year4,year4) = 1.000001*CAP.l(node,tec,year4,year4) ;
+        CAP.lo(node,tec,year4,year4) = 0.999999*CAP.l(node,tec,year4,year4) ;
+        ACT.fx(node,tec,year4,year4,mode,time) = ACT.l(node,tec,year4,year4,mode,time) ;
+        CAP_NEW_UP.fx(node,tec,year4) = CAP_NEW_UP.l(node,tec,year4) ;
+        CAP_NEW_LO.fx(node,tec,year4) = CAP_NEW_LO.l(node,tec,year4) ;
+        ACT_UP.fx(node,tec,year4,time) = ACT_UP.l(node,tec,year4,time) ;
+        ACT_LO.fx(node,tec,year4,time) = ACT_LO.l(node,tec,year4,time) ;
+
+
+        Display year,year4,year_all,year_all2,model_horizon ;
     ) ; # end of the recursive-dynamic loop
 
 ) ; # end of if statement for the selection betwen perfect-foresight or recursive-dynamic model
@@ -167,4 +182,4 @@ COST_NODAL_NET.L(node, year)$(NOT macro_base_period(year)) = (
             AND map_node(node2,node) AND cat_year(type_year,year) ),
         emission_scaling(type_emission,emission) * tax_emission(node2,type_emission,type_tec,type_year)
         * EMISS.L(node,emission,type_tec,year) )
-) ;
+) / 1000 ;

message_ix/model/MESSAGE/sets_maps_def.gms

@@ -152,6 +152,8 @@ Sets
     land_type       types of land use
     year_all        years (over entire model horizon)
     year (year_all) years included in a model instance (for myopic or rolling-horizon optimization)
+    year4 (year_all) years included in a model instance (for myopic or rolling-horizon optimization)
+    year_hist (year_all) historical year
     time            subannual time periods (seasons - days - hours)
     shares          share constraint relations
     relation        generic linear relations
@@ -301,6 +303,7 @@ Sets
     addon(tec)                              technologies that are an add-on to other (parent) technologies
     type_addon                              types of add-on technologies (that can be applied mutually exclusive)
     cat_addon(type_addon,addon)             mapping of add-on technologies to respective add-on technology types
+    newtec(tec)                             new technologies that will experience technological learning
     type_year                               types of year aggregations
     cat_year(type_year,year_all)            mapping of years to respective categories
     type_emission                           types of emission aggregations

message_ix/model/MESSAGE_master.gms

@@ -46,6 +46,7 @@ $ONGLOBAL
 
 ** scenario/case selection - this must match the name of the MsgData_<%%%>.gdx input data file **
 $SETGLOBAL data "<your datafile name here>"
+*$SETGLOBAL data "Westeros_Electrified_baseline"
 
 ** MACRO mode
 * "none": MESSAGEix is run in stand-alone mode
@@ -58,6 +59,12 @@ $SETGLOBAL macromode "none"
 * rolling horizon (period-by-period, recursive-dynamic with limited foresight - 'number of years of foresight'
 $SETGLOBAL foresight "0"
 
+** for recursive dynamic approach, this is to activate/deactivate technology learning module **
+* deactivate technology learning - 0
+* activate technology learning - 1
+$SETGLOBAL learningmode "0"
+
+
 ** add a comment and name extension for model report files (e.g. run-specific info, calibration notes) - optional **
 $SETGLOBAL comment ""