Scaling tool to improve LP matrix (cleaned)

message_ix/model/MESSAGE/data_load.gms

@@ -20,6 +20,7 @@ $LOAD balance_equality, time_relative
 $LOAD shares
 $LOAD addon, type_addon, cat_addon, map_tec_addon
 $LOAD storage_tec, level_storage, map_tec_storage
+$LOAD newtec
 
 * Version information; conditional load to allow older GDX files
 $ifthen gdxSetType ixmp_version
@@ -79,6 +80,13 @@ is_fixed_extraction, is_fixed_stock, is_fixed_new_capacity, is_fixed_capacity, i
 fixed_extraction, fixed_stock, fixed_new_capacity, fixed_capacity, fixed_activity, fixed_land
 * storage parameters
 storage_initial, storage_self_discharge, time_order
+
+* - - - - - - - - - - - - - -
+* learning related Parameters
+*** learning and economies of scale parameters
+alpha, beta_unit, beta_proj, gamma_unit, gamma_proj
+*** initial condition
+inv_cost_refidx, knref_unit, sizeref_unit, sizeref_proj
 ;
 
 

message_ix/model/MESSAGE/model_learningeos.gms

@@ -0,0 +1,171 @@
+
+PARAMETERS
+* inherrited parameters
+  cap_new2(newtec,year_all2)             'newly installed capacity of learning technologies'
+  bin_cap_new(newtec,year_all2)          'binary parameter of newly installed capacity'
+
+* learning and economies of scale parameters
+  alpha(newtec)                          'technology cost learning parameter                     n.a.'
+  beta_unit(newtec)                      'economy of scale parameter at unit level               n.a.'
+  beta_proj(newtec)                      'economy of scale parameter at project level            n.a.'
+  gamma_unit(newtec)                     'unit scale-up rate                                     GW per experience'
+  gamma_proj(newtec)                     'project scale-up rate                                  units per experience'
+
+* initial condition
+  inv_cost_refidx(newtec)                'initial capex (indexed to 1)                           n.a.'
+  knref_unit(newtec)                     'initial number of unit                                 units'
+  sizeref_unit(newtec)                   'initial size of unit                                   GW per unit'
+  sizeref_proj(newtec)                   'initial size of project                                units per project'
+
+* log2 parameters
+  log2_cap_new2(newtec,year_all2)        'log2 of newly installed capacity                       n.a.'
+  log2_inv_cost_refidx(newtec)           'log2 of indexed initial capex                          n.a.'
+  log2_knref_unit(newtec)                'log2 of initial number of unit                         n.a.'
+  log2_sizeref_unit(newtec)              'log2 of initial size of unit                           n.a.'
+  log2_sizeref_proj(newtec)              'log2 of initial size of project                        n.a.'
+;
+
+* calculating log2 parameters value
+log2_inv_cost_refidx(newtec)    = log2(inv_cost_refidx(newtec)) ;
+log2_knref_unit(newtec)         = log2(knref_unit(newtec)) ;
+log2_sizeref_unit(newtec)       = log2(sizeref_unit(newtec)) ;
+log2_sizeref_proj(newtec)       = log2(sizeref_proj(newtec)) ;
+
+* calculate the length of historical periods for indexing
+SCALAR hist_length                       'the length of historical periods' ;
+hist_length = card(year_all2) - card(model_horizon);
+
+VARIABLES
+* linear variables
+  IC(newtec,year_all2)                   'specific capital cost                                  dollar per kW'
+  N_UNIT(newtec,year_all2)               'number of new units                                    units'
+  KN_UNIT(newtec,year_all2)              'number of cumulative units                             units'
+  S_UNIT(newtec,year_all2)               'unit size                                              GW per unit or GtCO2 per unit'
+  S_PROJ(newtec,year_all2)               'project size                                           units per project'
+
+* log2 variables
+  LOG2_IC(newtec,year_all2)              'log2 of capital cost                                   n.a.'
+  LOG2_N_UNIT(newtec,year_all2)          'log2 of number of new units                            n.a.'
+  LOG2_KN_UNIT(newtec,year_all2)         'log2 of number of cumulative units                     n.a.'
+  LOG2_S_UNIT(newtec,year_all2)          'log2 of unit size                                      n.a.'
+  LOG2_S_PROJ(newtec,year_all2)          'log2 of project size                                   n.a.'
+
+* objective variable
+  OBJECT                                 'objective function'
+;
+
+* initializing lower bound of variable to avoid log2 singularity
+N_UNIT.lo(newtec,year_all2)  = 1E-15 ;
+KN_UNIT.lo(newtec,year_all2) = 1E-15 ;
+S_UNIT.lo(newtec,year_all2)  = sizeref_unit(newtec) ;
+S_PROJ.lo(newtec,year_all2)  = sizeref_proj(newtec) ;
+
+
+EQUATIONS
+  OBJECTIVE_INNER                        'objective value: total investment cost of learning technologies'
+  CAP_NEW_BALANCE                        'newly installed installed capacity balance'
+  KN_UNIT_LOG                            'log2 of cumulative number of units equality'
+  N_UNIT_LOG                             'log2 of number of units equality'
+  S_UNIT_LOG                             'log2 of unit size equality'
+  S_PROJ_LOG                             'log2 of project size equality'
+  CAPEX_ESTIMATE                         'estimating average capex'
+  CUMUL_UNIT_INI                         'calculate initial cumulative units'
+  CUMUL_UNIT                             'calculate cumulative units'
+  UNIT_SCALEUP_LIM_INI                   'calculate initial unit scale-up limit'
+  UNIT_SCALEUP_LIM                       'calculate unit scale-up limit'
+  PROJ_SCALEUP_LIM_INI                   'calculate initial project scale-up limit'
+  PROJ_SCALEUP_LIM                       'calculate project scale-up limit'
+  INV_COST_LOG                           'calculate investment cost'
+  UNIT_SIZELB                            'unit size lower bound'
+  UNIT_SIZE_NOBUILTYEAR                  'make unit size constant if no new capacity'
+  PROJ_SIZELB                            'project size lower bound'
+  PROJ_SIZE_NOBUILTYEAR                  'make project size constant if no new capacity'
+;
+
+
+OBJECTIVE_INNER..        OBJECT =e= sum((node,newtec,year_all2),
+                         IC(newtec,year_all2) * cap_new2(newtec,year_all2)) ;
+
+CAP_NEW_BALANCE(newtec,year_all2)$(bin_cap_new(newtec,year_all2) = 1)..
+         log2_cap_new2(newtec,year_all2) =e=
+         LOG2_N_UNIT(newtec,year_all2) + LOG2_S_UNIT(newtec,year_all2) ;
+
+KN_UNIT_LOG(newtec,year_all2)..
+         LOG2_KN_UNIT(newtec,year_all2) =e= log2(KN_UNIT(newtec,year_all2)) ;
+
+N_UNIT_LOG(newtec,year_all2)$(bin_cap_new(newtec,year_all2) = 1)..
+         LOG2_N_UNIT(newtec,year_all2) =e=
+         log2(N_UNIT(newtec,year_all2)) ;
+
+S_UNIT_LOG(newtec,year_all2)..
+         LOG2_S_UNIT(newtec,year_all2) =e= log2(S_UNIT(newtec,year_all2)) ;
+
+S_PROJ_LOG(newtec,year_all2)..
+         LOG2_S_PROJ(newtec,year_all2) =e= log2(S_PROJ(newtec,year_all2)) ;
+
+CAPEX_ESTIMATE(newtec,year_all2)..  LOG2_IC(newtec,year_all2) =e=
+         LOG2_IC(newtec,year_all2-1)
+         - alpha(newtec) * [LOG2_KN_UNIT(newtec,year_all2) - LOG2_KN_UNIT(newtec,year_all2-1)] * bin_cap_new(newtec,year_all2)
+         - beta_unit(newtec) * [LOG2_S_UNIT(newtec,year_all2) - LOG2_S_UNIT(newtec,year_all2-1)] * bin_cap_new(newtec,year_all2)
+         - beta_proj(newtec) * [LOG2_S_PROJ(newtec,year_all2) - LOG2_S_PROJ(newtec,year_all2-1)] * bin_cap_new(newtec,year_all2) ;
+
+CUMUL_UNIT_INI(newtec,year_all2)$(ord(year_all2) le (hist_length+1))..
+         KN_UNIT(newtec,year_all2) =e=
+         knref_unit(newtec) + N_UNIT(newtec,year_all2) * bin_cap_new(newtec,year_all2) ;
+
+CUMUL_UNIT(newtec,year_all2)$(ord(year_all2) gt (hist_length+1))..
+         KN_UNIT(newtec,year_all2) =e=
+         KN_UNIT(newtec,year_all2-1) + N_UNIT(newtec,year_all2) * bin_cap_new(newtec,year_all2) ;
+
+UNIT_SCALEUP_LIM_INI(newtec,year_all2)$(ord(year_all2) le (hist_length+1))..
+         S_UNIT(newtec,year_all2) =l=
+         sizeref_unit(newtec)
+         + gamma_unit(newtec) * [LOG2_KN_UNIT(newtec,year_all2) - log2_knref_unit(newtec)] ;
+
+UNIT_SCALEUP_LIM(newtec,year_all2)$(ord(year_all2) gt (hist_length+1))..
+         S_UNIT(newtec,year_all2) =l=
+         S_UNIT(newtec,year_all2-1)
+         + gamma_unit(newtec) * [LOG2_KN_UNIT(newtec,year_all2) - LOG2_KN_UNIT(newtec,year_all2-1)] ;
+
+PROJ_SCALEUP_LIM_INI(newtec,year_all2)$(ord(year_all2) le (hist_length+1))..
+         S_PROJ(newtec,year_all2) =l=
+         sizeref_proj(newtec)
+         + gamma_proj(newtec) * [LOG2_KN_UNIT(newtec,year_all2) - log2_knref_unit(newtec)] ;
+
+PROJ_SCALEUP_LIM(newtec,year_all2)$(ord(year_all2) gt (hist_length+1))..
+         S_PROJ(newtec,year_all2) =l=
+         S_PROJ(newtec,year_all2-1)
+         + gamma_proj(newtec) * [LOG2_KN_UNIT(newtec,year_all2) - LOG2_KN_UNIT(newtec,year_all2-1)] ;
+
+INV_COST_LOG(newtec,year_all2)..
+         IC(newtec,year_all2) =e= 2**LOG2_IC(newtec,year_all2) ;
+
+UNIT_SIZELB(newtec,year_all2)..
+         S_UNIT(newtec,year_all2) =g= sizeref_unit(newtec) ;
+
+UNIT_SIZE_NOBUILTYEAR(newtec,year_all2)$(bin_cap_new(newtec,year_all2) = 0 and
+         ord(year_all2) gt (hist_length+1))..
+         S_UNIT(newtec,year_all2) =e= S_UNIT(newtec,year_all2-1) ;
+
+PROJ_SIZELB(newtec,year_all2)..
+         S_PROJ(newtec,year_all2) =g= sizeref_proj(newtec) ;
+
+PROJ_SIZE_NOBUILTYEAR(newtec,year_all2)$(bin_cap_new(newtec,year_all2) = 0 and
+         ord(year_all2) gt (hist_length+1))..
+         S_PROJ(newtec,year_all2) =e= S_PROJ(newtec,year_all2-1) ;
+
+
+* declaring model equations
+* please keep model equations listed in this format
+* 'all' statement will include all MESSAGE-ix equations and increase learning module solution time
+model learningeos /
+  OBJECTIVE_INNER,  CAP_NEW_BALANCE,  KN_UNIT_LOG,  N_UNIT_LOG,  S_UNIT_LOG,
+  S_PROJ_LOG, CAPEX_ESTIMATE, CUMUL_UNIT_INI, CUMUL_UNIT,  UNIT_SCALEUP_LIM_INI,
+  UNIT_SCALEUP_LIM, PROJ_SCALEUP_LIM_INI, PROJ_SCALEUP_LIM, INV_COST_LOG,
+  UNIT_SIZELB,
+  UNIT_SIZE_NOBUILTYEAR,
+  PROJ_SIZELB,
+  PROJ_SIZE_NOBUILTYEAR
+  /;
+* keep this option for testing purpose
+*option nlp = minos;

message_ix/model/MESSAGE/model_setup.gms

@@ -21,6 +21,14 @@ $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"
 
+** include scaler commands
+$IF NOT SET scaler $SETGLOBAL scaler "MsgScaler_Default"
+
+** 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 +92,5 @@ $INCLUDE MESSAGE/scaling_investment_costs.gms
 *----------------------------------------------------------------------------------------------------------------------*
 
 $INCLUDE MESSAGE/model_core.gms
+$INCLUDE MESSAGE/model_learningeos.gms
+$INCLUDE scaler/%scaler%.gms

message_ix/model/MESSAGE/model_solve.gms

@@ -6,6 +6,18 @@
 * including the required accounting of investment costs beyond the model horizon.
 ***
 
+Parameter
+    count_iter                           'iteration counter'
+    prev_OBJ                             'previous objective value'
+    delta_OBJ                            'difference between current and previous objective value'
+    inv_cost_ini(node,newtec)            'initial investment cost assumption at first model year'
+;
+
+count_iter = 1;
+prev_OBJ   = 0;
+delta_OBJ  = 1;
+inv_cost_ini(node,newtec) = sum(first_period, inv_cost(node,newtec,first_period));
+
 if (%foresight% = 0,
 ***
 * Perfect-foresight model
@@ -22,6 +34,45 @@ if (%foresight% = 0,
 * include all model periods in the optimization horizon (excluding historical periods prior to 'first_period')
     year(year_all)$( model_horizon(year_all) ) = yes ;
 
+    IF(%learningmode% = 1,
+        put_utility 'log' /'+++ Solve the perfect-foresight with learning version of MESSAGEix +++ ' ;
+        option threads = 4 ;
+        while(count_iter <= 10 and
+              delta_OBJ >= 0.01,
+              Solve MESSAGE_LP using LP minimizing OBJ ;
+
+*             passing CAP_NEW values to update cap_new2 data for unit and size optimization,
+*             make bin parameters equal to 1 when technology is built, and 0 if otherwise,
+*             and passing log2_cap_new parameter to learning module
+*             NOTE: there is no log 0, thus binary parameter is added as helper to avoid error.
+*                   multiplication with the binary parameter in learning module negates effect of this helper
+              cap_new2(newtec,year_all2) = sum(node, CAP_NEW.l(node,newtec,year_all2));
+              bin_cap_new(newtec,year_all2) = sum(node, CAP_NEW.l(node,newtec,year_all2));
+              bin_cap_new(newtec,year_all2)$(bin_cap_new(newtec,year_all2) > 1E-10) = 1 ;
+              log2_cap_new2(newtec,year_all2) = log2(sum(node, CAP_NEW.l(node,newtec,year_all2)) + [1-bin_cap_new(newtec,year_all2)] ) ;
+              Solve learningeos using nlp minimizing OBJECT;
+
+*             update inv_cost values using indexed (normalized) cost IC
+              inv_cost(node,newtec,year_all2) = IC.l(newtec,year_all2) * inv_cost_ini(node,newtec);
+              if(count_iter = 1,
+                      delta_OBJ = 1 ;
+              else
+                      delta_OBJ = (prev_OBJ - OBJ.l)/prev_OBJ ;
+              );
+              prev_OBJ = OBJ.l ;
+              display count_iter, delta_OBJ;
+              display year,year_all,year_all2,model_horizon,inv_cost_ini ;
+              count_iter = count_iter + 1 ;
+        );
+
+    ELSE
+*       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 ;
+        display year,year_all,year_all2,model_horizon,inv_cost_ini ;
+    );
+
 * write a status update to the log file, solve the model
     put_utility 'log' /'+++ Solve the perfect-foresight version of MESSAGEix +++ ' ;
     Solve MESSAGE_LP using LP minimizing OBJ ;
@@ -96,6 +147,7 @@ else
 * 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 ;
+        year4(year_all2)$((ord(year_all2) < ord(year_all))) = 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 '  +++ ' ;
@@ -115,7 +167,55 @@ else
             ABORT "MESSAGEix did not solve to optimality!"
         ) ;
 
+        IF(%learningmode% = 1,
+*       passing CAP_NEW values to update cap_new2 data for unit and size optimization,
+*       make bin parameters equal to 1 when technology is built, and 0 if otherwise,
+*       and passing log2_cap_new parameter to learning module
+*       NOTE: there is no log 0, thus binary parameter is added as helper to avoid error.
+*       multiplication with the binary parameter in learning module negates effect of this helper
+        cap_new2(newtec,year_all2) = sum(node, CAP_NEW.l(node,newtec,year_all2));
+        bin_cap_new(newtec,year_all2) = sum(node, CAP_NEW.l(node,newtec,year_all2));
+        bin_cap_new(newtec,year_all2)$(bin_cap_new(newtec,year_all2) > 1E-10) = 1 ;
+        log2_cap_new2(newtec,year_all2) = log2(sum(node, CAP_NEW.l(node,newtec,year_all2)) + [1-bin_cap_new(newtec,year_all2)] ) ;
+
+        Solve learningeos using nlp minimizing OBJECT;
+
+*       update inv_cost values using indexed (normalized) cost IC
+        inv_cost(node,newtec,year_all2) = IC.l(newtec,year_all2) * inv_cost_ini(node,newtec);
+
+        display hist_length, bin_cap_new, IC.l, inv_cost, cap_new2, CAP_NEW.l,inv_cost_ini;
+        );
+
 * fix all variables of the current iteration period 'year_all' to the optimal levels
+        ext_fix(node,commodity,grade,year4) =  EXT.l(node,commodity,grade,year4) ;
+        cap_new_fix(node,newtec,year4) = CAP_NEW.l(node,newtec,year4) ;
+        cap_fix(node,newtec,year4,year4) = CAP.l(node,newtec,year4,year4) ;
+        act_fix(node,newtec,year4,year4,mode,time) = ACT.l(node,newtec,year4,year4,mode,time) ;
+
+        ext_fix(node,commodity,grade,year4)$(ext_fix(node,commodity,grade,year4) > 1E-10) =  1 ;
+        cap_new_fix(node,newtec,year4)$(cap_new_fix(node,newtec,year4) > 1E-10) = 1 ;
+        cap_fix(node,newtec,year4,year4)$(cap_fix(node,newtec,year4,year4) > 1E-10) = 1 ;
+        act_fix(node,newtec,year4,year4,mode,time)$(act_fix(node,newtec,year4,year4,mode,time) > 1E-10) = 1 ;
+
+*        EXT.fx(node,commodity,grade,year4) =  EXT.l(node,commodity,grade,year4) ;
+        CAP_NEW.fx(node,newtec,year4)$(cap_new_fix(node,newtec,year4) = 1) = CAP_NEW.l(node,newtec,year4) ;
+*        CAP.fx(node,newtec,year4,year4) = CAP.l(node,newtec,year4,year4) ;
+        ACT.fx(node,newtec,year4,year4,mode,time)$(act_fix(node,newtec,year4,year4,mode,time) = 1) = ACT.l(node,newtec,year4,year4,mode,time) ;
+*        CAP_NEW_UP.fx(node,newtec,year4) = CAP_NEW_UP.l(node,newtec,year4) ;
+*        CAP_NEW_LO.fx(node,newtec,year4) = CAP_NEW_LO.l(node,newtec,year4) ;
+*        ACT_UP.fx(node,newtec,year4,time) = ACT_UP.l(node,newtec,year4,time) ;
+*        ACT_LO.fx(node,newtec,year4,time) = ACT_LO.l(node,newtec,year4,time) ;
+* relaxed "fix"
+*        EXT.up(node,commodity,grade,year4) = 1.000001*EXT.l(node,commodity,grade,year4) ;
+*        EXT.lo(node,commodity,grade,year4) = 0.999999*EXT.l(node,commodity,grade,year4) ;
+**        CAP_NEW.up(node,newtec,year4) = 1.000001*CAP_NEW.l(node,newtec,year4) ;
+**        CAP_NEW.lo(node,newtec,year4) = 0.999999*CAP_NEW.l(node,newtec,year4) ;
+*        CAP.up(node,newtec,year4,year4) = 1.000001*CAP.l(node,newtec,year4,year4) ;
+*        CAP.lo(node,newtec,year4,year4) = 0.999999*CAP.l(node,newtec,year4,year4) ;
+**        ACT.up(node,newtec,year4,year4,mode,time) = 1.000001*ACT.l(node,newtec,year4,year4,mode,time) ;
+**        ACT.lo(node,newtec,year4,year4,mode,time) = 0.999999*ACT.l(node,newtec,year4,year4,mode,time) ;
+
+$ontext
         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) ;
@@ -125,7 +225,9 @@ else
         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) ;
-
+$offtext
+
+        Display year,year4,year_all,year_all2,model_horizon,inv_cost_ini ;
     ) ; # end of the recursive-dynamic loop
 
 ) ; # end of if statement for the selection betwen perfect-foresight or recursive-dynamic model