(!******************************************************
Mosel R Example Problems
========================
file energyforecast.mos
```````````````````````
Forecast of energy demand using regression analysis
on historic data
(c) 2015 Fair Isaac Corporation
author: S.Lannez, Jul. 2015
*******************************************************!)
model forecast
options noimplicit
uses "mmsheet" ! Spreadsheet I/O
uses "mmsystem"
uses "r" ! Load R interface
parameters
DATADIR=""
! Historical data
HISTFILE=DATADIR+'energyforecast.dat'
! Install R packages on the-fly
R_SETUP = false
! The peak threshold. If the ratio between
! the daily demand and the previous year average
! daily demand is greater than the following threshold
! then the period is considered a peak period.
FORECAST_PEAK_THR = 2.5
! The peak period probability threshold. If
! the predicted probability of being in a
! peak period is greater than the following
! threshold, then the period is considered
! to be a peak period.
FORECAST_PEAK_PROB = 0.5
end-parameters
! **** Historical data **** !
declarations
Days: set of string ! Days
Years: set of integer ! Years
Periods: range ! Time periods
HourTypes: set of integer ! Hour types
DayTypes: set of integer ! Day types
WeekTypes: set of integer ! Week types
AttrGDP: set of string ! YearGDP attributes
AttrIncome: set of string ! Income attributes
HourDemand: dynamic array(Days,Periods) of real ! Historical demand
YearGDP: dynamic array(Years,AttrGDP) of real ! Historical YearGDP
YearIncome: dynamic array(Years,AttrIncome) of real ! Historical income
first,last: date
tnow: datetime
end-declarations
! **** Return the hour type ****
function gethourtype(t: integer): integer
declarations
h=t/48 ! Input data has 48 time periods
end-declarations
returned := floor(h*HourTypes.size)
end-function
! **** Return the week type ****
function getweektype(d: date): integer
returned := getdaynum(date(d)) div 7
end-function
! **** Return the day type ****
function getdaytype(d: date): integer
declarations
daynum: integer
Summer=91..274
SummerHoliday=268..274
WinterHoliday=84..90
end-declarations
returned := 1 ! unknown
daynum := getdaynum(d)
if daynum in Summer then ! Summer
if daynum in SummerHoliday then
returned := 2 ! Summer public holiday
else
case getweekday(d) of
1,2,3,4,5: returned := 0 ! Summer week
6,7: returned := 1 ! Summer weekend
end-case
end-if
else ! Winter
if daynum in WinterHoliday then
returned := 5 ! Winter public holiday
else
case getweekday(d) of
1,2,3,4,5: returned := 3 ! Winter week
6,7: returned := 4 ! Winter weekend
end-case
end-if
end-if
end-function
! **** Run a regression analysis of forecast the energy consumption ****
procedure forecast(first,last: date,predict:boolean)
! Sets the parameters of the forecast
! function
declarations
TimePeriods,_TimePeriods: set of string
! Derived data
PeriodAvgEnergy: dynamic array(_TimePeriods) of real ! Average daily energy of the year
! Observed data or forecasted from external tool
PeriodGDP: dynamic array(TimePeriods) of real ! Current year GDP
PeriodIncome: dynamic array(TimePeriods) of real ! Current year income
PeriodWeekType: dynamic array(TimePeriods,WeekTypes) of boolean ! Day type
PeriodDayType: dynamic array(TimePeriods,DayTypes) of boolean ! Day type
PeriodHourType: dynamic array(TimePeriods,HourTypes) of boolean ! Hour type
PeriodWeekTypeS: dynamic array(TimePeriods) of string ! String representation of week type
PeriodDayTypeS: dynamic array(TimePeriods) of string ! String representation of day type
PeriodHourTypeS: dynamic array(TimePeriods) of string ! String representation of hour type
! Observed demand
PeriodDemand: dynamic array(TimePeriods) of real ! baseline demand
PeriodPeak: dynamic array(TimePeriods) of real ! peak demand
PeriodIPeak: dynamic array(TimePeriods) of real ! peak period
! Period forecast
PeriodDemandP: dynamic array(TimePeriods) of real ! Baseline energy demand forecast
PeriodPeakP: dynamic array(TimePeriods) of real ! Peak energy demand forecast
PeriodIPeakP: dynamic array(TimePeriods) of real ! Probability of being in a peak
! Misc
id,newdata: string
i,y: integer
v,w: real
E,_E,__E,err,a,b: real
Calendar: dynamic array(Years,Days,Periods) of string
end-declarations
writeln('-'*60)
if (predict) then
writeln('Forecasting...')
else
writeln('Running regression...')
if (R_SETUP) then
Reval('install.packages("mlogit",repos="http://cran.r-project.org")')
end-if
end-if
! Create lookup array
writeln("Creating calendar...")
forall(d in Days, p in Periods | exists(HourDemand(d,p))) do
Calendar(getyear(date(d)),d,p) := d+"-"+if(p>=10,""+p,"0"+p)
end-do
! Compute previous year average consumed energy.
! The forecast function is based on energy demand growth
! and needs relative demand increase from year-to-year
forall(yy in Years) do
i := 0
E := sum(d in Days, p in Periods | exists(Calendar(yy-1,d,p)),i as counter) HourDemand(d,p)
if (i>0) then
E := E/i
else
E := 0
end-if
forall(d in Days, p in Periods | exists(Calendar(yy,d,p))) do
PeriodAvgEnergy(Calendar(yy,d,p)) := E
end-do
end-do
! Populate the predictors used in the regression or
! during the prediction.
writeln("Populating predictors...")
forall(d in Days | date(d)>=first and date(d)<=last) do
y := getyear(date(d))
forall (p in Periods | exists(Calendar(y,d,p))) do
id := Calendar(y,d,p)
! Continous input
PeriodGDP(id) := YearGDP(y,'YearGDPGrowthRate')
PeriodIncome(id) := YearIncome(y,'Income')
! Continuous input
PeriodWeekType(id,getweektype(date(d))) := true
PeriodDayType(id,getdaytype(date(d))) := true
PeriodHourType(id,gethourtype(p)) := true
! Discrete input
PeriodWeekTypeS(id) := ""+getweektype(date(p))
PeriodDayTypeS(id) := ""+getdaytype(date(p))
PeriodHourTypeS(id) := ""+gethourtype(p)
! Value to be forecasted
if (exists(HourDemand(d,p))) then
w := HourDemand(d,p) / PeriodAvgEnergy(id)
! We split baseline and peak forecasting
! in two different models in order
! to improve the quality of the prediction
if (w>=FORECAST_PEAK_THR) then
PeriodPeak(id) := w-FORECAST_PEAK_THR
PeriodIPeak(id) := 1
else
PeriodPeak(id) := 0
PeriodIPeak(id) := 0
end-if
PeriodDemand(id) := w
end-if
end-do
end-do
writeln("Loading data to R...")
! Copy data to R
Rset('timeperiods',TimePeriods)
Reval('MyData <- data.frame(row.names=timeperiods)')
Rset('MyData$gdp',PeriodGDP)
Rset('MyData$income',PeriodIncome)
Rset('MyData$weektype',PeriodWeekType)
Rset('MyData$daytype',PeriodDayType)
Rset('MyData$hourtype',PeriodHourType)
Rset('MyData$weektypestr',PeriodWeekTypeS)
Rset('MyData$daytypestr',PeriodDayTypeS)
Rset('MyData$hourtypestr',PeriodHourTypeS)
! Fill NAs cell with 0
Reval('MyData$hourtype[ is.na(MyData$hourtype) ] <- 0')
Reval('MyData$daytype[ is.na(MyData$daytype) ] <- 0')
Reval('MyData$weektype[ is.na(MyData$weektype) ] <- 0')
! When not predicting we need to set the demand, peak
! and indicator peak period in order to run the regression
if (not predict) then
Rset('MyData$demand',PeriodDemand)
Rset('MyData$peak',PeriodPeak)
Rset('MyData$ipeak',PeriodIPeak)
end-if
! Predict
if (not predict) then
! Define our prediction model
! This model is more advanced and takes into account the income and
! gdp of the year, the weektype, daytype and hourtype. It will generates
! a set of coefficients that can be used to forecast the demand by taking
! into account the forecasted mean temperature, gdp and annual income of
! the future periods.
writeln("Fitting '"+TimePeriods.size+"' time periods...")
! Fitting the peak
writeln("... peak indicators ...")
Reval('peakModelI <- glm(ipeak ~ weektypestr + daytypestr + hourtypestr, data=MyData, family="binomial")')
! Uncomment to get significance of each parameter.
! This computation can take minutes.
! Rprint('anova(peakModelI,test="Chisq")') ! Print predictors relevance
! Fitting the baseline consumption
! We assume baseline is a linear response of the gdp/income/peak probability
writeln("... baseline consumption ...")
Reval('demandModel <- glm(demand ~ gdp + income + weektype + daytype*hourtype, data=MyData)')
! Uncomment to get significance of each parameter.
! This computation can take minutes.
! Rprint('anova(demandModel,test="Chisq")') ! Print predictors relevance
! Fitting the peak consumption
! We assume peak consumption is an exponential response of the gdp/daytype/hourtype predictors
writeln("... peak consumption ...")
Reval('peakModel <- glm(peak ~ gdp + income + weektype + daytype*hourtype, data=MyData)')
! Uncomment to get significance of each parameter.
! This computation can take minutes.
! Rprint('anova(peakModel,test="Chisq")') ! Print predictors relevance
! Save the forecasting model
! Uncomment the following line to extract the models to
! a file that can be read by R
! Reval('save(demandModel, peakModel, peakModelI, file="'+DATADIR+'forecaster.Rdata")') ! Debug
! Free some memory
Reval('rm(MyData)')
else
writeln("Predicting '"+TimePeriods.size+"' time periods...")
! Predicting
Reval('MyData$ipeak <- predict(peakModelI, newdata = MyData, type = "response")') ! Peak hours
Reval('MyData$demand <- predict(demandModel, newdata = MyData)') ! Baseline
Reval('MyData$peak <- predict(peakModel, newdata = MyData)') ! Peak value
writeln("==== Summary ===========")
writeln("= Peak period forecast =")
Rprint('summary(MyData$ipeak)')
writeln("= Baseline demand forecast =")
Rprint('summary(MyData$demand)')
writeln("= Peak demand forecast =")
Rprint('summary(MyData$peak)')
! Copy the data to Mosel
Rgetarr('mmarray(MyData,"demand")',PeriodDemandP)
Rgetarr('mmarray(MyData,"ipeak")',PeriodIPeakP)
Rgetarr('mmarray(MyData,"peak")',PeriodPeakP)
! Clean up R environment
Reval('rm(MyData)')
writeln("==== Analysis ==============")
! Print some summary
E := 0 ; _E := 0; __E := 0 ; err := 0
forall(t in TimePeriods) do
__E += PeriodAvgEnergy(t)
_E += PeriodDemand(t) * PeriodAvgEnergy(t)
E += PeriodDemandP(t) * PeriodAvgEnergy(t)
err += abs(PeriodDemandP(t) - PeriodDemand(t))*PeriodAvgEnergy(t)
end-do
writeln(" === Baseline forecaster === ")
writeln("Previous year total energy : ",textfmt(__E/10e6,6,3)," TWh")
writeln("Original total energy : ",textfmt(_E/10e6,6,3)," TWh")
writeln("Forecasted total energy : ",textfmt(E/10e6,6,3)," TWh")
if (maxlist(E,_E)>0) then
writeln("Absolute error : ",textfmt(round(abs(E-_E)/10e3),6)," GWh")
writeln("Relative error : ",textfmt(round(100*abs(E-_E)/maxlist(E,_E)),6)+" %")
end-if
writeln("Cumulated period error : ",textfmt(err/10e3,6)," GWh")
E := 0 ; _E := 0; __E := 0 ; err := 0
forall(t in TimePeriods) do
__E += PeriodAvgEnergy(t)
_E += PeriodPeak(t) * PeriodAvgEnergy(t)
E += PeriodPeakP(t) * PeriodAvgEnergy(t)
err += abs(PeriodPeakP(t) - PeriodPeak(t))*PeriodAvgEnergy(t)
end-do
writeln(" === Peak forecaster === ")
writeln("Previous year total energy : ",textfmt(__E/10e6,6,3)," TWh")
writeln("Original total energy : ",textfmt(_E/10e6,6,3)," TWh")
writeln("Forecasted total energy : ",textfmt(E/10e6,6,3)," TWh")
if (maxlist(E,_E)>0) then
writeln("Absolute error : ",textfmt(round(abs(E-_E)/10e3),6)," GWh")
writeln("Relative error : ",textfmt(round(100*abs(E-_E)/maxlist(E,_E)),6)+" %")
end-if
writeln("Cumulated period error : ",textfmt(err/10e3,6)," GWh")
end-if
writeln("Done.")
end-procedure
! ************************ Application ************************ !
! Define the sets
AttrGDP := {'YearGDP','YearGDPGrowthRate'} ! Available YearGDP data
AttrIncome := {'Income','Population','PopulationGrowthRate','IncomeGrowthRate'} ! Available income data
Periods := 0..47 ! Set of time periods in a day
HourTypes := 0..47 ! Set of hour types
DayTypes := {-1,0,1,2,3,4,5} ! Set of day types
! Helper tool to convert from data.frame to
! simple arrays
Reval('mmarray <- function(ar,col) { a <- cbind(ar[,col]) ; colnames(a) <- col ; rownames(a) <- rownames(ar) ; a[,col] }')
! In this example we have historical data
! for the years 2010,2011,2012,2013. We have
! gdp and income forecast for 2014 and we
! want to derive the energy demand.
!
! Note: today is the 1st of January, 2014
!
! Load historical and forecast data.
! All data before 1st of January, 2014 are considered
! past data
initializations from HISTFILE
HourDemand as "E_DEMAND" ! HourDemand
YearGDP as "IND_GDP" ! YearGDP
YearIncome as "IND_INCOME" ! Income
end-initializations
! **** Run regression (2011-2013) **** !
first := date('2011-01-01')
last := date('2013-12-31')
forecast(first,last,false)
! **** Predict demand (2014) **** !
first := date('2014-02-01')
last := date('2014-04-30')
forecast(first,last,true)
end-model
|
