Many changes, see below

Backbone: added sets and parameters for energy modelling for buildings
Plots: Added option -o (output file name) to do_r_plot.r
Using a alternative objective file

.gitignore

@@ -29,3 +29,12 @@ EfficiencyPiecewise.xlsx
 to_UCD/
 *.zip
 sandbox/
+.Rproj.user
+/.Rhistory
+/gamsproject.dat
+/howToRead.txt
+/inc/1c_parameters.gms.bak
+/inc/1e_inputs.gms.bak
+/plots/.RData
+/plots/.Rhistory
+*.Rproj

backbone_definitions_4plots.gms

+* Define empty sets, parameters and variables (copy pasted from backbone)
+
+* Check current GAMS version
+$ife %system.gamsversion%<240 $abort GAMS distribution 24.0 or later required!
+
+* Set default debugging level
+$if not set debug $setglobal debug 0
+
+* Default values for input and output dir as well as input data GDX file and index sheet when importing data from Excel file
+* When reading an Excel file, you can opt to read the file only if the Gdxxrw detects changes by using 'checkDate' for
+*   input_excel_checkdate. It is off by default, since there has been some problems with it.
+$if not set input_dir $setglobal input_dir 'input'
+$if not set output_dir $setglobal output_dir 'output'
+$if not set input_file_gdx $setglobal input_file_gdx 'inputData.gdx'
+$if not set input_excel_index $setglobal input_excel_index 'INDEX'
+$if not set input_excel_checkdate $setglobal input_excel_checkdate ''
+
+* Make sure output dir exists
+$if not dexist '%output_dir%' $call 'mkdir %output_dir%'
+
+* Activate end of line comments and set comment character to '//'
+$oneolcom
+$eolcom //
+$onempty   // Allow empty data definitions
+
+* Output file streams
+Files log /''/, gdx /''/, f_info /'%output_dir%/info.txt'/;
+
+* Include options file to control the solver (if it does not exist, uses defaults)
+$ifthen exist '%input_dir%/1_options.gms'
+    $$include '%input_dir%/1_options.gms';
+$endif
+
+* === Libraries ===============================================================
+$libinclude scenred2
+
+* === Definitions, sets, parameters and input data=============================
+$include 'inc/1a_definitions.gms'   // Definitions for possible model settings
+$include 'inc/1b_sets.gms'          // Set definitions used by the models
+$include 'inc/1c_parameters.gms'    // Parameter definitions used by the models
+$include 'inc/1d_results.gms'       // Parameter definitions for model results
+$include 'inc/1e_inputs.gms'        // Load input data
\ No newline at end of file

inc/1b_sets.gms

@@ -179,6 +179,14 @@ $if defined scenario
     gnGroup(grid, node, group) "Combination of grids and nodes in particular groups"
     sGroup(s, group) "Samples in particular groups"
 
+* --- Energy modelling for buildings ---------------------------------------------------------------
+    node_building
+    node_building_DHWT(node),
+    node_building_envelope_mass(node),
+    node_building_interior_air_and_furniture(node),
+    node_building_internal_mass(node)
+    node2node_building(node_building,node)
+
 * --- Set of timeseries that will be read from files between solves -----------
     mTimeseries_loop_read(mType, timeseries) "Those time series that will be read between solves"
 ;

inc/1c_parameters.gms

@@ -211,3 +211,9 @@ Parameters
     p_stepLengthNoReset(mType, f, t) "Length of an interval in hours - includes also lengths of previously realized intervals"
     p_s_discountFactor(s) "Discount factor for samples when using a multi-year horizon"
 ;
+
+* --- Energy modelling for buildings -----------------------------------------
+Parameters
+    ts_priceElspot(t) "Elspot prices (EUR/MWh)"
+    building_squares(node_building) "Building square meters (m2)" 
+;

inc/1e_inputs.gms

@@ -90,6 +90,11 @@ $loaddc p_groupPolicyUnit
 $loaddc p_groupPolicyEmission
 $loaddc gnss_bound
 $loaddc uss_bound
+* Energy modelling for buildings
+$loaddc ts_priceElspot
+$loaddc node_building node_building_DHWT node_building_envelope_mass node_building_interior_air_and_furniture node_building_internal_mass
+$loaddc node2node_building
+$loaddc building_squares
 
 $gdxin
 

input/1_options.gms

@@ -2,12 +2,17 @@ Options
 
 // Solution gap: the first one reached will end iteration
 optca = 0       // Absolute gap between the found solution and the best possible solution
-optcr = 0.005  // Relative gap between the found solution and the best possible solution
+optcr = 0.00  // Relative gap between the found solution and the best possible solution
 
 // Solver options
 Solvelink = %Solvelink.ChainScript%
 
 threads = 1          // How many cores the solver can use: 0 = all cores; negative values = all cores - n
-solprint = silent   // Don't print a lot of stuff into .lst?
+*solprint = silent   // Don't print a lot of stuff into .lst?
+option limrow=9, limcol=99;
+option solver=cbc;
 ;
 
+$onInclude
+$onSymXRef;
+

input/2c_alternative_objective.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
+
+
+* =============================================================================
+* --- Objective Function Definition -------------------------------------------
+* =============================================================================
+Parameter
+  ts_priceElspotNode(node,t);
+ts_priceElspotNode(node,t)=1E6;
+*node_building_DHWT node_building_envelope_mass node_building_interior_air_and_furniture node_building_internal_mass
+*ts_priceElspotNode(node_building_DHWT,t)=ts_priceElspot(t);
+*ts_priceElspotNode(node_building_interior_air_and_furniture,t)=ts_priceElspot(t);
+ts_priceElspotNode("link_node",t)=ts_priceElspot(t);
+
+
+q_obj ..
+
+    + v_obj * 1e6 // Objective function valued in MEUR instead of EUR (or whatever monetary unit the data is in)
+
+    =E=
+
+    // Sum over all the samples, forecasts, and time steps in the current model
+    + sum(msft(m, s, f, t),
+        // Probability (weight coefficient) of (s,f,t)
+        + p_msft_probability(m, s, f, t)
+            * p_s_discountFactor(s) // Discount costs
+            * [
+                // Time step length dependent costs
+                + p_stepLength(m, f, t)
+                    * [
+                        // Dummy variable penalties
+                        // Energy balance feasibility dummy varible penalties
+                        +   sum(inc_dec,
+                            + sum(gn(grid, node),
+                                + vq_gen(inc_dec, grid, node, s, f, t)
+                                    * ts_priceElspotNode(node,t)
+* original                          *( PENALTY_BALANCE(grid, node)${not p_gnBoundaryPropertiesForStates(grid, node, 'balancePenalty', 'useTimeSeries')}
+* original                          + ts_node_(grid, node, 'balancePenalty', s, f, t)${p_gnBoundaryPropertiesForStates(grid, node, 'balancePenalty', 'useTimeSeries')}
+* original                                    )
+                                ) // END sum(gn)
+                            ) // END sum(inc_dec)
+                        ] // END * p_stepLength
+                ]  // END * p_sft_probability(s,f,t)
+        ) // END sum over msft(m, s, f, t)
+;

input/inputDataAdjustments.gms

@@ -40,34 +40,61 @@ $gdxin
 Scalars
   modify_input      "Modify input (0=no, 1=yes: to set new limits)"                     / 1      /
   c2k               "Convert celcius to kelvin degrees"                                 / 275.15 /
-  start_temp_air    "Start value and lower limit for inside air temperature in Celcius" / 20     /
-  temp_up_air       "Upper limit for inside air temperature in Celcius"                 / 24     /
 ;
 Sets
- node_DHWT(node) /"IDA_ESBO_DH1.DHWT1000", "IDA_ESBO_DH2.DHWT1000", "IDA_ESBO_AB.DHWT1000"/
- node_envelope_mass(node) / "IDA_ESBO_DH1.envelope_mass", "IDA_ESBO_DH2.envelope_mass", "IDA_ESBO_AB.envelope_mass" /
- node_interior_air_and_furniture(node) / "IDA_ESBO_DH1.interior_air_and_furniture", "IDA_ESBO_DH2.interior_air_and_furniture", "IDA_ESBO_AB.interior_air_and_furniture" /
- node_internal_mass(node) / "IDA_ESBO_DH1.internal_mass", "IDA_ESBO_DH2.internal_mass", "IDA_ESBO_AB.internal_mass" /
+ dir   "Direction: lo ,reference ,up"           / lo, ref, up/
+ place "Measurement place in building"          / DHWT, envelope_mass, internal_mass, interior_air_and_furniture /
+ node_building                                  / IDA_ESBO_DH1, IDA_ESBO_DH2, IDA_ESBO_AB  /
+ node_building_DHWT(node)                       / "IDA_ESBO_DH1.DHWT1000", "IDA_ESBO_DH2.DHWT1000", "IDA_ESBO_AB.DHWT1000" /
+ node_building_envelope_mass(node)              / "IDA_ESBO_DH1.envelope_mass", "IDA_ESBO_DH2.envelope_mass", "IDA_ESBO_AB.envelope_mass" /
+ node_building_internal_mass(node)              / "IDA_ESBO_DH1.internal_mass", "IDA_ESBO_DH2.internal_mass", "IDA_ESBO_AB.internal_mass" /
+ node_building_interior_air_and_furniture(node) / "IDA_ESBO_DH1.interior_air_and_furniture", "IDA_ESBO_DH2.interior_air_and_furniture", "IDA_ESBO_AB.interior_air_and_furniture" /
+ node2node_building(node_building,node) /
+   IDA_ESBO_DH1.("IDA_ESBO_DH1.DHWT1000", "IDA_ESBO_DH1.envelope_mass", "IDA_ESBO_DH1.internal_mass","IDA_ESBO_DH1.interior_air_and_furniture")
+   IDA_ESBO_DH2.("IDA_ESBO_DH2.DHWT1000", "IDA_ESBO_DH2.envelope_mass", "IDA_ESBO_DH2.internal_mass","IDA_ESBO_DH2.interior_air_and_furniture")
+   IDA_ESBO_AB.("IDA_ESBO_AB.DHWT1000", "IDA_ESBO_AB.envelope_mass", "IDA_ESBO_AB.internal_mass","IDA_ESBO_AB.interior_air_and_furniture")
+ /;
+
+
+
+
+Table
+  temp(place,dir)   "Limit temperatures in Celcius"
+                             lo      ref        up
+DHWT                         60      60         80
+envelope_mass                0       20         50
+interior_air_and_furniture   20      20         24
+internal_mass                20      20         24
 ;
 
+*TODO
+Parameter building_squares (node_building) Building square meters /
+"IDA_ESBO_DH1"  135.56
+"IDA_ESBO_DH2"  145.33
+"IDA_ESBO_AB"  1608.19
+/;
+
 if(modify_input=1,
 * Defintion: p_gnBoundaryPropertiesForStates(grid,node,param_gnBoundaryTypes,param_gnBoundaryProperties)
-  display grid,node,node_DHWT,node_envelope_mass,node_interior_air_and_furniture,node_internal_mass,param_gnBoundaryTypes,param_gnBoundaryProperties,p_gnBoundaryPropertiesForStates;
-* Domestic Hot Water Tank (DHWT) is set to lower limit value at start
-  p_gnBoundaryPropertiesForStates('building',node_DHWT,'reference','useConstant') = 1;
-  p_gnBoundaryPropertiesForStates('building',node_DHWT,'reference','constant') = p_gnBoundaryPropertiesForStates('building',node_DHWT,'downwardLimit','constant') ;
-* set lower limit of temperature to 20C
-  p_gnBoundaryPropertiesForStates('building',node_envelope_mass,'downwardLimit','constant') = 275.15+20;
-  p_gnBoundaryPropertiesForStates('building',node_interior_air_and_furniture,'downwardLimit','constant') = 275.15+20;
-  p_gnBoundaryPropertiesForStates('building',node_internal_mass,'downwardLimit','constant') = 275.15+20;
-* set start value (reference) to lower limit
-  p_gnBoundaryPropertiesForStates('building',node_envelope_mass,'reference','constant') = p_gnBoundaryPropertiesForStates('building',node_envelope_mass,'downwardLimit','constant') ;
-  p_gnBoundaryPropertiesForStates('building',node_interior_air_and_furniture,'reference','constant') = p_gnBoundaryPropertiesForStates('building',node_interior_air_and_furniture,'downwardLimit','constant');
-  p_gnBoundaryPropertiesForStates('building',node_internal_mass,'reference','constant') = p_gnBoundaryPropertiesForStates('building',node_internal_mass,'downwardLimit','constant');
-* set upper limit for building temperature
-  p_gnBoundaryPropertiesForStates('building',node_envelope_mass,'upwardLimit','constant') = 275.15+20;
-  p_gnBoundaryPropertiesForStates('building',node_interior_air_and_furniture,'upwardLimit','constant') = 275.15+20;
-  p_gnBoundaryPropertiesForStates('building',node_internal_mass,'upwardLimit','constant') = 275.15+20;
+  display grid,node,node_building_DHWT,node_building_envelope_mass,node_building_interior_air_and_furniture,node_building_internal_mass,param_gnBoundaryTypes,param_gnBoundaryProperties,p_gnBoundaryPropertiesForStates;
+* Set temperature limits for Domestic Hot Water Tank (DHWT)
+  p_gn(grid, node_building_DHWT, 'boundStart') = 1;
+  p_gnBoundaryPropertiesForStates('building',node_building_DHWT,'reference','useConstant') = 1;
+  p_gnBoundaryPropertiesForStates('building',node_building_DHWT,'downwardLimit','constant') = c2k + temp('DHWT' ,'lo');
+  p_gnBoundaryPropertiesForStates('building',node_building_DHWT,'reference','constant') = c2k + temp('DHWT' ,'ref');
+  p_gnBoundaryPropertiesForStates('building',node_building_DHWT,'upwardLimit','constant')= c2k + temp('DHWT' ,'up');
+* Set temperature limits for building_envelope_mass in Celcius
+  p_gnBoundaryPropertiesForStates('building',node_building_envelope_mass,'downwardLimit','constant') = c2k + temp('envelope_mass' ,'lo ');
+  p_gnBoundaryPropertiesForStates('building',node_building_envelope_mass,'reference','constant') = c2k + temp('envelope_mass' ,'ref');
+  p_gnBoundaryPropertiesForStates('building',node_building_envelope_mass,'upwardLimit','constant') = c2k +temp('envelope_mass' ,'up');;
+* Set temperature (Celcius) limits for building_internal_mass in Celcius
+  p_gnBoundaryPropertiesForStates('building',node_building_internal_mass,'downwardLimit','constant') = c2k + temp('internal_mass','lo');
+  p_gnBoundaryPropertiesForStates('building',node_building_internal_mass,'reference','constant') = c2k + temp('internal_mass','ref');
+  p_gnBoundaryPropertiesForStates('building',node_building_internal_mass,'upwardLimit','constant') = c2k + temp('internal_mass','up');
+* Set temperature (Celcius) limits for building_interior_air_and_furniture in Celcius
+  p_gnBoundaryPropertiesForStates('building',node_building_interior_air_and_furniture,'downwardLimit','constant') = c2k + temp('interior_air_and_furniture' ,'lo');
+  p_gnBoundaryPropertiesForStates('building',node_building_interior_air_and_furniture,'reference','constant') = c2k + temp('interior_air_and_furniture' ,'ref');
+  p_gnBoundaryPropertiesForStates('building',node_building_interior_air_and_furniture,'upwardLimit','constant') = c2k + temp('interior_air_and_furniture' ,'up');
 ) ;
 
 
@@ -83,7 +110,10 @@ $gdxIn input\elspot_prices_2013.gdx
 $load tsIso elspotIsoBB ts_priceElspot=elspotBB
 $gdxin
 
+* Adding exogenous commodity electricity that can be bought
+
+
 * Write input with domains at execution time to take in consideration modifications
 execute_unloaddi "input\inputDataAdjusted1.gdx";
-execute_unloaddi "input\inputDataAdjusted1Squuezed.gdx" ts_priceElspot ts_influx effLevelGroupUnit p_gnBoundaryPropertiesForStates p_gnu_io p_gn p_unit node p_gnn unit unitUnittype grid p_s_discountFactor unittype
+execute_unloaddi "input\inputDataAdjusted1Squuezed.gdx" node2node_building ts_priceElspot ts_influx effLevelGroupUnit p_gnBoundaryPropertiesForStates p_gnu_io p_gn p_unit node p_gnn unit unitUnittype grid p_s_discountFactor unittype
 

output/info.txt

@@ -2,17 +2,17 @@
  MODEL RUN DETAILS                                                   
 
 User                tltoni
-Date                12/28/21
-Time                12:02:16
+Date                01/07/22
+Time                14:16:58
 GAMS version        37.1
 GAMS system         GAMS Base Module 37.1.0 r07954d5 Released Nov 11, 2021 WEI x86 64bit/MS Window
-version             v1.5-58-g4ec63ce+
+version             v1.5-60-g04e5406+
 
              time (s)
 
-Compilation         1
-Execution           8
-Total              32
+Compilation         0
+Execution          19
+Total              52
 
 
  MODEL FEATURES                                                      

r/do_r_plot.gms → plots/do_r_plot.gms

 $title Run Rscript do_r_plot.r
 Set
- t "Time (e.g. hours)       " /t1*t200/
+ t "Time (e.g. hours)       "               /t1*t200/
  c "Structure measurement point (Celcius) " / inside, wall, outside/
 ;
 Parameter
 p(t,c) "Building temperature (celcius) "
 ;
 p(t,c)=ord(t)*ord(t)*ord(c);
-execute_unload "%gams.curDir%r\plot4r.gdx";
-execute 'Rscript "%gams.curDir%r\do_r_plot.r" -h'
-execute 'Rscript "%gams.curDir%r\do_r_plot.r" -f "%gams.curDir%r\plot4r.gdx" -p p -w 50 -a TRUE';
+execute_unload "%gams.wDir%plots\plot4r.gdx";
+execute 'Rscript "%gams.wDir%plots\do_r_plot.r" -h'
+execute 'Rscript "%gams.wDir%plots\do_r_plot.r" -f "%gams.wDir%plots\plot4r.gdx" -p p -w 50 -a TRUE';
 scalar myerrorlevel;
 myerrorlevel = errorlevel;
 if(myerrorlevel=0,
-  execute.async 'SumatraPDF.exe  "%gams.curDir%Rplots.pdf"';
+  execute.async 'SumatraPDF.exe  "%gams.wDir%plots\Rplots.pdf"';
 )

r/do_r_plot.r → plots/do_r_plot.r

@@ -24,41 +24,47 @@ library(tools)
 library(rlang)
 library(forcats)
 
-#Preprocessing : delete old pdf
-pdf_name="Rplots.pdf"
-tryCatch({
-  if(file.exists(pdf_name)){
-    invisible(file.remove(pdf_name))
-  }
-}
-  , error = function(ex) { print(ex) ; stop("FAILED: see above message", call.=FALSE) }  
-)
 
 #CMD Line Option Management and Error Checking
 option_list = list(
-  make_option(c("-f", "--file"), type="character", default="plot4r.gdx", 
-              help="GDX file name [default= %default]", metavar="character"),
+  make_option(c("-i", "--input"), type="character", default="plot4r.gdx", 
+              help="Input file name (GDX) [default= %default]", metavar="character"),
+  make_option(c("-o", "--output"), type="character", default="Rplots.pdf", 
+              help="Output file name (PDF)  [default= %default]", metavar="character"),
   make_option(c("-p", "--parameter"), type="character", default="p", 
               help="parameter name in gdx, expects 3-dimensional, e.g. p(x,y) [default= %default]", metavar="character"),
   make_option(c("-w", "--window"), type="integer", default="", 
               help="do additional plots by splitting x (x-axis) to the specified time step length [default= %default]", metavar="integer"),
   make_option(c("-a", "--archive"), type="logical", default=FALSE, 
-              help="archive the created pdf in folder ggplots [default= %default]", metavar="logical")
+              help="archive the created pdf in folder archive [default= %default]", metavar="logical")
     ); 
 opt_parser = OptionParser(option_list=option_list, description = program_description)
 opt = parse_args(opt_parser);
 
+#Preprocessing : set working directory and delete old pdf
+r_dir = "/Users/tltoni/GIT/FlexiB/backbone_building_v2/plots"
+setwd(r_dir)
+pdf_name=opt$output
+tryCatch({
+  if(file.exists(pdf_name)){
+    invisible(file.remove(pdf_name))
+  }
+}
+, error = function(ex) { print(ex) ; stop("FAILED: see above message", call.=FALSE) }  
+)
+
+
 # check that gdx file and gdx reader is found
-if (is.null(opt$file) || isFALSE(file.exists(opt$file))) {
+if (is.null(opt$input) || isFALSE(file.exists(opt$input))) {
   print_help(opt_parser)
-  stop(sprintf("Input GDX file not found: %s/%s",getwd(),opt$file), call.=FALSE)
+  stop(sprintf("Input GDX file not found: %s/%s",getwd(),opt$input), call.=FALSE)
 }
 if (! require(gdxrrw))      stop ("gdxrrw package is not available")
 if (0 == igdx(silent=TRUE)) stop ("the gdx shared library has not been loaded")
 
 # read gdx parameter p of p(x,y)
 tryCatch({
-  gdxName <- opt$file
+  gdxName <- opt$input
   symName <- opt$parameter
   df_p  <- rgdx.param(gdxName,symName, ts=TRUE)
   }
@@ -74,7 +80,7 @@ if(ncol(df_p)!=3){
 # read set x and y from parameter p(x,y)
 xvar_desc <- vector("character")
 tryCatch({
-  gdxName <- opt$file
+  gdxName <- opt$input
   symName <- opt$parameter
   xvar=as_string(sym(colnames(df_p)[1]))
   df_x <- rgdx.set(gdxName, xvar, names=xvar, ts=TRUE)
@@ -84,7 +90,7 @@ tryCatch({
 )
 yvar_desc <- vector("character")
 tryCatch({
-  gdxName <- opt$file
+  gdxName <- opt$input
   symName <- opt$parameter
   yvar=as_string(sym(colnames(df_p)[2]))
   df_y <- rgdx.set(gdxName, yvar, names=yvar, ts=TRUE)
@@ -154,7 +160,7 @@ if(file.exists(pdf_name)){
   print(sprintf("%s",pdf_path))
   flush.console()
   if(opt$archive){
-    dir_name="ggplots"
+    dir_name="archive"
     if(!dir.exists(dir_name)){
       dir.create(dir_name)
     }

plots/plot_ts_influx.gms

+$title Plot backbone output with R
+$ontext
+Plot backbone output with R
+
+1. Reads backbones full output, e.g. debug.gdx
+2. Uses do_r_plot.r do create figures of seleceted parameters
+
+Toni Lastusilta (VTT) , 2021-12-31
+$offtext
+
+Set
+  t      "Model time steps (in hours)"
+  node   "Track exogenous commodities (MWh/h)"
+  nodeP  "Track exogenous commodities by m2 (MWh/h)"
+  node_building
+  node2node_building(node_building,node)
+;
+
+Parameter
+  ext(t,node)   "External power (weather) influencing building"
+  ext2(t,nodeP) "External power (weather) influencing building - normalized by bulding m2"
+;
+
+$batinclude backbone_definitions_4plots.gms
+$onlisting
+
+* CREATE PLOTS
+$onmultiR
+$gdxin output\archive\out_alt_obj.gdx
+$load node nodeP=node node_building
+$load ts_influx  node2node_building building_squares
+$gdxin
+ext(t,node)=sum((grid,f)$(ord(f)=1 and sameas(grid,"building")),ts_influx(grid,node,f,t)$(ord(f)=1 and ord(grid)=1));
+ext2(t,nodeP)=sum((node2node_building(node_building,node))$sameas(node,nodeP),ext(t,node)/building_squares(node_building));
+execute_unload 'plots/plot4r.gdx', t, node, nodeP, ext, ext2;
+execute 'Rscript "%gams.wDir%plots\do_r_plot.r" -i "%gams.wDir%plots/plot4r.gdx" -o "%gams.wDir%output/ts_influx1.pdf"  -p ext -w 1000 -a FALSE';
+scalar myerrorlevel;
+myerrorlevel = errorlevel;
+if(myerrorlevel=0,
+  execute.async 'SumatraPDF.exe  "%gams.wDir%output/ts_influx1.pdf"';
+)
+
+execute 'Rscript "%gams.wDir%plots\do_r_plot.r" -i "%gams.wDir%plots/plot4r.gdx" -o "%gams.wDir%output/ts_influx2.pdf"  -p ext2 -w 1000 -a FALSE';
+myerrorlevel = errorlevel;
+if(myerrorlevel=0,
+  execute.async 'SumatraPDF.exe  "%gams.wDir%output/ts_influx2.pdf"';
+)
+
+
+
+
+
+
+
+
+
+
+
+
+
+

plots/simple_r_plot_create_data.gms

+*Create data for simple Rscript create_simple_r_plot.r
+set
+ i "Set i " /i1*i5/
+ j "Set j " /j1*j3/
+;
+parameter
+p(i,j) "parameter p"
+;
+p(i,j)=ord(i)*ord(j);
+execute_unload "%gams.wDir%\plots\simple_r_plot_data.gdx";
+*execute 'Rscript "%gams.curDir%plots\create_simple_r_plot.r"';
+

plots/simple_r_plot_create_plot.r

+
+# A simple example to plot a GDX file with R
+library(reshape2)
+library(plyr)
+library(gdxrrw)
+library(ggplot2)
+
+r_dir = "/Users/tltoni/GIT/FlexiB/backbone_building_v2/plots"
+# we use gams backbone project directory to avoid path issues
+gams_project_dir = "/Users/tltoni/GIT/FlexiB/backbone_building_v2"
+gdx_filename = "simple_r_plot_data.gdx"
+setwd(r_dir)
+if (! require(gdxrrw))      stop ("gdxrrw package is not available")
+if (0 == igdx(silent=TRUE)) stop ("the gdx shared library has not been loaded")
+
+if(!file.exists(gdx_filename)){
+  exec_string= sprintf("gams plots/simple_r_plot_create_data.gms wdir=%s",gams_project_dir)
+  system(exec_string)
+}
+
+#gdxInfo("gdx_filename")
+tryCatch({
+  gdxName <- gdx_filename
+  symName <- "p"
+  df <- rgdx.param(gdxName,symName)
+#  ix <- rgdx.set(gdxName,'i',names="I")
+#  jx <- rgdx.set(gdxName,'j',names="J")
+}
+, error = function(ex) { stop(ex) ; FALSE }
+)
+
+g  <- ggplot(df, aes(x=i, y=p, group=j, color=j)) + geom_line() 
+plot(g)
+