{ "type": "Article", "authors": [ { "type": "Person", "familyNames": [ "Madsen" ], "givenNames": [ "MikkelL.Sørensen,JanK.MøllerandHenrik" ] } ], "description": [ { "type": "Paragraph", "content": [ "For electricity grid operators, the planning of grid operations depends on having accurate models for forecasting wind\npower production on a number of different time resolutions. Together, these time resolutions create what has become\nknown as a temporal hierarchy. Previous studies have considered methods of reconciling forecast hierarchies inspired by\nleast squares methodologies that produce coherent and more accurate forecasts. In this study, we highlight some\nchallenges in the established approach when applied to wind power production, and consider methods which more\nappropriately take into account the full conditional probability densities. We suggest methods using maximum likelihood\ntechniques to estimate the prediction variance ahead of time. Using base forecasts from a commercial forecast provider\ntogether with simpler forecasting models, we test the modified approach against the established reconciliation approach\non data from Danish wind farms. The results show significant improvements in accuracy when compared to both the\nstate-of-the-art commercial forecasts and the simpler models." ] } ], "identifiers": [], "references": [ { "type": "Article", "id": "bib-bib1", "authors": [], "title": "\nG. Athanasopoulos, R. A. Ahmed and R. J. Hyndman,\nHierarchical\nforecasts for Australian domestic tourism. Int. J. Forecast.\n25, 146–166 (2009) ", "url": "https://dx.doi.org/10.1016/j.ijforecast.2008.07.004" }, { "type": "Article", "id": "bib-bib2", "authors": [], "title": "\nG. Athanasopoulos, R. J. Hyndman, N. Kourentzes and F. Petropoulos,\nForecasting with\ntemporal hierarchies. Eur. J. Oper. Res. 262, 60–74\n(2017) ", "url": "https://dx.doi.org/10.1016/j.ejor.2017.02.046" }, { "type": "Article", "id": "bib-bib3", "authors": [], "title": "\nH. G. Bergsteinsson, J. K. Møller, P. Nystrup, Ó. P. Pálsson, D. Guericke and\nH. Madsen, Heat load\nforecasting using adaptive temporal hierarchies. Appl. Energ.\n292, article no. 116872 (2021) ", "url": "https://dx.doi.org/10.1016/j.apenergy.2021.116872" }, { "type": "Article", "id": "bib-bib4", "authors": [], "title": "\nF. Cribari-Neto and A. Zeileis,\nBeta regression in R.\nJ. Stat. Software 34, 1–24 (2010) ", "url": "https://dx.doi.org/10.18637/jss.v034.i02" }, { "type": "Article", "id": "bib-bib5", "authors": [], "title": "\nC. Di Modica, P. Pinson and S. Ben Taieb,\nOnline forecast\nreconciliation in wind power prediction. Electr. Power Syst. Res.\n190, article no. 106637 (2021) ", "url": "https://dx.doi.org/10.1016/j.epsr.2020.106637" }, { "type": "Article", "id": "bib-bib6", "authors": [], "title": "\nM. E. Hansen, P. Nystrup, J. K. Møller and H. Madsen, Reconciliation of wind\npower forecasts in spatial hierarchies. Wind Energ. 26,\n615–632 (2023) ", "url": "https://doi.org/10.1002/we.2819" }, { "type": "Article", "id": "bib-bib7", "authors": [], "title": "\nR. J. Hyndman, R. A. Ahmed, G. Athanasopoulos and H. L. Shang,\nOptimal combination\nforecasts for hierarchical time series. Comput. Statist. Data Anal.\n55, 2579–2589 (2011) ", "url": "https://dx.doi.org/10.1016/j.csda.2011.03.006" }, { "type": "Article", "id": "bib-bib8", "authors": [], "title": "\nR. J. Hyndman,\nA. J. Lee and E. Wang,\nFast computation of\nreconciled forecasts for hierarchical and grouped time series.\nComput. Statist. Data Anal.\n97, 16–32 (2016) ", "url": "https://dx.doi.org/10.1016/j.csda.2015.11.007" }, { "type": "Article", "id": "bib-bib9", "authors": [], "title": "\nN. Kourentzes and G. Athanasopoulos,\nCross-temporal\ncoherent forecasts for Australian tourism. Ann. Tourism\nRes. 75, 393–409 (2019) ", "url": "https://dx.doi.org/10.1016/j.annals.2019.02.001" }, { "type": "Article", "id": "bib-bib10", "authors": [], "title": "\nM. Lange, On the uncertainty of\nwind power predictions – Analysis of the forecast accuracy and statistical\ndistribution of errors. J. Sol. Energ. Eng. 127, 177–184 (2005) ", "url": "https://dx.doi.org/10.1115/1.1862266" }, { "type": "Article", "id": "bib-bib11", "authors": [], "title": "\nH. Madsen and P. Thyregod, Introduction to general and generalized linear\nmodels. CRC Press, Boca Raton, FL (2011) ", "url": "https://doi.org/10.1201/9781439891148" }, { "type": "Article", "id": "bib-bib12", "authors": [], "title": "\nJ. K. Møller, M. Zugno and H. Madsen,\nProbabilistic forecasts of wind\npower generation by stochastic differential equation models. J. Forecast.\n35, 189–205 (2016) ", "url": "https://dx.doi.org/10.1002/for.2367" }, { "type": "Article", "id": "bib-bib13", "authors": [], "title": "\nH. A. Nielsen, H. Madsen and T. S. Nielsen,\nUsing quantile regression to extend\nan existing wind power forecasting system with probabilistic forecasts.\nWind Energ. 9, 95–108 (2006) ", "url": "https://dx.doi.org/10.1002/we.180" }, { "type": "Article", "id": "bib-bib14", "authors": [], "title": "\nP. Nystrup, E. Lindström, P. Pinson and H. Madsen,\nTemporal hierarchies\nwith autocorrelation for load forecasting. Eur. J. Oper. Res.\n280, 876–888 (2020) ", "url": "https://dx.doi.org/10.1016/j.ejor.2019.07.061" }, { "type": "Article", "id": "bib-bib15", "authors": [], "title": "\nP. Pinson, Wind energy: Forecasting\nchallenges for its operational management. Stat. Sci. 28,\n564–585 (2013) ", "url": "https://dx.doi.org/10.1214/13-STS445" }, { "type": "Article", "id": "bib-bib16", "authors": [], "title": "\nM. L. Sørensen, P. Nystrup, M. B. Bjerregård, J. K. Møller, P. Bacher\nand H. Madsen, Recent developments in multivariate wind and solar power\nforecasting. WIREs Energy Environ.\n12, article no. e465 (2023) ", "url": "https://doi.org/10.1002/wene.465" }, { "type": "Article", "id": "bib-bib17", "authors": [], "title": "\nS. L. Wickramasuriya, G. Athanasopoulos and R. J. Hyndman,\nOptimal forecast\nreconciliation for hierarchical and grouped time series through trace\nminimization. J. Am. Stat. Assoc. 114, 804–819 (2019)\n", "url": "https://dx.doi.org/10.1080/01621459.2018.1448825" }, { "type": "Article", "id": "bib-bib18", "authors": [], "title": "\nS. L. Wickramasuriya, B. A. Turlach and R. J. Hyndman,\nOptimal non-negative\nforecast reconciliation. Stat. Comput. 30, 1167–1182 (2020)\n", "url": "https://dx.doi.org/10.1007/s11222-020-09930-0" } ], "title": "Reconciling temporal hierarchies of wind power production with forecast-dependent variance structures", "meta": {}, "content": [ { "type": "Heading", "id": "Sx1", "depth": 1, "content": [ "Introduction" ] }, { "type": "Paragraph", "id": "Sx1.p1", "content": [ "As part of the ongoing green transition, the share of renewable energy being incorporated into the electricity grid is\nincreasing. With this comes the challenge of renewable energy sources, such as wind power, being stochastic in nature,\ni.e., power production depends on, e.g., the wind speed. Thus, in order to plan and control the electricity grid,\ntransmission system operators (TSOs) require reliable and accurate forecasts of power production. With targets such as\nthe European Union having 70% of power production stemming from renewable energy sources\nby 2050,", { "type": "Note", "id": "idm11", "noteType": "Footnote", "content": [ { "type": "Paragraph", "id": "footnote1", "content": [ "IRENA, Global energy transformation: A roadmap to 2050;\n", { "type": "Link", "target": "https://www.irena.org/publications/2019/Apr/Global-energy-transformation-A-roadmap-to-2050-2019Edition", "content": [ "https://www.irena.org/publications/2019/Apr/Global-energy-transformation-A-roadmap-to-2050-2019Edition" ] }, "." ] } ] }, "\nthis challenge is only becoming greater. Hence, the importance of accurate and reliable wind power\nforecasts is more crucial now than ever before." ] }, { "type": "Paragraph", "id": "Sx1.p2", "content": [ "The main goal in this respect is thus to reduce the uncertainty of the forecasted power production as much as possible.\nOne of the more recent developments within forecasting, which has been shown to reduce the forecast uncertainty, is the\nmethod of hierarchical forecast reconciliation. The method in its current form was suggested by [", { "type": "Cite", "target": "bib-bib7", "content": [ "7" ] }, "],\nwith further developments by [", { "type": "Cite", "target": "bib-bib8", "content": [ "8" ] }, ", ", { "type": "Cite", "target": "bib-bib17", "content": [ "17" ] }, ", ", { "type": "Cite", "target": "bib-bib14", "content": [ "14" ] }, "].\nHierarchical forecast reconciliation utilizes the hierarchical structures that are often naturally present in data by ensuring that\nforecasts on different resolutions are coherent. Coherency refers to the fact that the layers should add up throughout\nthe hierarchy in a natural way. An example of this comes from the study on Australian domestic tourism by [", { "type": "Cite", "target": "bib-bib1", "content": [ "1" ] }, "],\nwhere forecasts are produced at the national scale, state scale, zone scale, and regional\nscale. These forecasts form a hierarchy, where one would expect that all the regional forecasts within a zone add up to\nthe zone forecast, all the zone forecasts within a state add up to the state forecast, etc. When forecasts add up in\nthis manner, the forecasts are said to be coherent. However, when the forecasts on each layer are produced\nindependently of each other this will not be the case. The purpose of hierarchical forecast reconciliation is thus\nto produce a set of reconciled forecasts that have the coherency property based on a set of incoherent base forecasts." ] }, { "type": "Paragraph", "id": "Sx1.p3", "content": [ "Hierarchies can be either spatial, as was the case for the example on Australian domestic tourism, temporal, as has\nbeen explored by [", { "type": "Cite", "target": "bib-bib2", "content": [ "2" ] }, "] and [", { "type": "Cite", "target": "bib-bib14", "content": [ "14" ] }, "], or spatio-temporal, e.g., [", { "type": "Cite", "target": "bib-bib9", "content": [ "9" ] }, "].\nIn a spatial hierarchy one might examine regional forecasts, see, e.g., [", { "type": "Cite", "target": "bib-bib7", "content": [ "7" ] }, "]. In a temporal hierarchy,\nwhich will be the case for this study, the forecasts should add up based on the timescale such that, e.g., the three\nmonthly forecasts within a quarter add up to quarterly forecasts, and the four quarterly forecasts add up to yearly\nforecasts." ] }, { "type": "Paragraph", "id": "Sx1.p4", "content": [ "Forecast reconciliation of temporal hierarchies has shown great promise in reducing the uncertainty of the base\nforecasts and producing more accurate forecasts at all layers [", { "type": "Cite", "target": "bib-bib14", "content": [ "14" ] }, "]. However, the majority of studies on\ntemporal hierarchies only show such improvements using simple models to produce the base forecasts, such as the ARIMA\nand exponential smoothing models by [", { "type": "Cite", "target": "bib-bib2", "content": [ "2" ] }, "]. These models also often require constant data\navailability, which under real-time operational conditions cannot always be guaranteed. The established reconciliation\napproach also assumes that the covariance structure between the forecasts is constant in time, which for many\napplications may not be the case. For wind power production it is well-known that the forecast uncertainty is highly\ndependent on, e.g., the wind speed [", { "type": "Cite", "target": "bib-bib13", "content": [ "13" ] }, "]. Furthermore, wind power production is naturally bounded between\nzero and the installed capacity, which poses a challenge when reconciling, since the established method of reconciling\nassumes unbounded data. These challenges call for more refined methods to generate reconciled forecasts." ] }, { "type": "Paragraph", "id": "Sx1.p5", "content": [ "In the literature, forecast reconciliation has been applied previously to cases of both wind power forecasting and\nsolar power forecasting, which both face similar issues. The paper [", { "type": "Cite", "target": "bib-bib16", "content": [ "16" ] }, "] provides a brief overview\nof the literature on this. In [", { "type": "Cite", "target": "bib-bib6", "content": [ "6" ] }, "] the authors test some of the current methods for forecast\nreconciliation in a spatial hierarchy with a similar dataset to our case study and are able to achieve significant\nimprovements in accuracy. They do not, however, account for the issue of the non-constant variance structure. The need\nfor varying and mean-dependent variance when reconciling forecasts has also been highlighted by [", { "type": "Cite", "target": "bib-bib3", "content": [ "3" ] }, "]\nfor the case of heat load forecasts, although the approach therein mostly addresses slower\nseasonal variations. Additionally, [", { "type": "Cite", "target": "bib-bib18", "content": [ "18" ] }, "] demonstrates a method of changing the forecast\nreconciliation methodology to account for non-negativity." ] }, { "type": "Paragraph", "content": [ "The purpose of the present study is thus:" ] }, { "type": "Paragraph", "content": [] }, { "type": "List", "items": [ { "type": "ListItem", "content": [ { "type": "Paragraph", "id": "Sx1.I1.i1.p1", "content": [ "To use real measurements and weather forecasts to examine whether improvements to commercial state-of-the-art\nwind power production forecasts can be achieved through forecast reconciliation." ] } ] }, { "type": "ListItem", "content": [ { "type": "Paragraph", "id": "Sx1.I1.i2.p1", "content": [ "To propose adjustments to the established reconciliation approach (MinT-Shrink by [", { "type": "Cite", "target": "bib-bib17", "content": [ "17" ] }, "]),\nwhich addresses the level and time-dependent covariance structure in the wind power framework.\n" ] } ] } ], "order": "Unordered", "meta": { "listType": "bullet" } }, { "type": "Paragraph", "id": "Sx1.p9", "content": [ "This study is structured such that Section ", { "type": "Cite", "target": "S1", "content": [ "1" ] }, " introduces the concept of forecast reconciliation.\nSection ", { "type": "Cite", "target": "S2", "content": [ "2" ] }, " presents the wind power production data used in the study. Section ", { "type": "Cite", "target": "S3", "content": [ "3" ] }, " then\nexamines how base forecasts are produced, and introduces the proposed adjustments to the reconciliation process.\nSection ", { "type": "Cite", "target": "S5", "content": [ "5" ] }, " shows the results, and lastly Section ", { "type": "Cite", "target": "S6", "content": [ "6" ] }, " concludes on the findings and\ngives some further perspectives." ] }, { "type": "Heading", "id": "S1", "depth": 1, "content": [ "1 Forecast reconciliation" ] }, { "type": "Paragraph", "id": "S1.p1", "content": [ "For the temporal hierarchy in this study, prediction horizons up to 24 hours will be examined. So, to introduce the\nconcepts of reconciling forecasts based on temporal hierarchies, a simple example hierarchy is constructed within the\n24-hour prediction frame with 6-hour, 12-hour and 24-hour forecast resolutions. This example could, e.g., cover the\ntotal wind power production in the stated intervals.\nThe hierarchical structure of this example is illustrated in Figure ", { "type": "Cite", "target": "S1-F1", "content": [ "1" ] }, "." ] }, { "type": "Figure", "id": "S1-F1", "caption": [ { "type": "Paragraph", "content": [ "Temporal hierarchy with three levels covering a 24-hour period." ] } ], "content": [ { "type": "ImageObject", "contentUrl": "172-f1.png", "mediaType": "image/png", "meta": { "inline": false } } ] }, { "type": "Paragraph", "content": [ "As illustrated by the figure, for such a hierarchy of forecasts to be coherent, the bottom-level 6-hour forecasts\nshould add up pairwise to the second level 12-hour forecast, and the two 12-hour forecasts should add up to the highest\nlevel daily forecast. This summation property of the hierarchy is described by a summation matrix ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝐒∈ℝn×m", "meta": { "altText": "\\mathbf{S}\\in\\mathbb{R}^{n\\times m}" } }, ", where ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "n", "meta": { "altText": "n" } }, " is the number of individual base forecasts (in this case 7) and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "m", "meta": { "altText": "m" } }, " is the number\nof forecasts on the lowest resolution (in this case 4). For the hierarchy in Figure ", { "type": "Cite", "target": "S1-F1", "content": [ "1" ] }, ", the summation\nmatrix is\ngiven by:" ] }, { "type": "MathBlock", "id": "S1.Ex1", "mathLanguage": "mathml", "text": "𝐒=[1111110000111000010000100001].", "meta": { "altText": "\\mathbf{S}=\\begin{bmatrix}1&1&1&1\\\\\n1&1&0&0\\\\\n0&0&1&1\\\\\n1&0&0&0\\\\\n0&1&0&0\\\\\n0&0&1&0\\\\\n0&0&0&1\\end{bmatrix}." } }, { "type": "Paragraph", "content": [ "The rows of the summation matrix define how the four bottom-level forecasts add up to each level of the hierarchy, so,\ne.g., the top-level forecast should be equal to the sum of all the four bottom-level forecasts. This also implies that\nthe base forecasts (", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y^∈ℝn", "meta": { "altText": "\\hat{y}\\in\\mathbb{R}^{n}" } }, ") are ordered from top\nto bottom, i.e." ] }, { "type": "MathBlock", "id": "S1.E1", "label": "(1)", "mathLanguage": "mathml", "text": "y^=[y^24⁢h24⁢hy^24⁢h12⁢hy^12⁢h12⁢hy^24⁢h6⁢hy^18⁢h6⁢hy^12⁢h6⁢hy^6⁢h6⁢h].\\notag", "meta": { "altText": "\\hat{y}=\\begin{bmatrix}\\hat{y}^{24h}_{24h}\\\\[2.84526pt]\n\\hat{y}^{12h}_{24h}\\\\[2.84526pt]\n\\hat{y}^{12h}_{12h}\\\\[2.84526pt]\n\\hat{y}^{6h}_{24h}\\\\[2.84526pt]\n\\hat{y}^{6h}_{18h}\\\\[2.84526pt]\n\\hat{y}^{6h}_{12h}\\\\[2.84526pt]\n\\hat{y}^{6h}_{6h}\\end{bmatrix}.\\notag" } }, { "type": "Paragraph", "content": [ "For clarity, the notation in this study uses superscripts for naming or other distinctions and subscripts for indices,\nso, e.g., ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y^18⁢h6⁢h", "meta": { "altText": "\\hat{y}^{6h}_{18h}" } }, " is the forecast on the 6-hour resolution 13 to 18 hours ahead." ] }, { "type": "Paragraph", "id": "S1.p3", "content": [ "With the summation matrix and the base forecasts as defined above, the most commonly used method of forecast\nreconciliation, the MinT-Shrink [", { "type": "Cite", "target": "bib-bib17", "content": [ "17" ] }, "], can be introduced." ] }, { "type": "Paragraph", "content": [ "The method is based on a generalized least squares approach, which minimizes the coherency errors, i.e., the difference\nbetween the reconciled forecast (", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y~∈ℝn", "meta": { "altText": "\\tilde{y}\\in\\mathbb{R}^{n}" } }, ") and the base forecast (", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y^", "meta": { "altText": "\\hat{y}" } }, "). This implies solving\nthe following optimization\nproblem" ] }, { "type": "MathBlock", "id": "S1.Ex2", "mathLanguage": "mathml", "text": "miny~(y~−y^)T⁢𝚺⁢(y~−y^)subject toy~=𝐒𝐆⁢y~,", "meta": { "altText": "\\min_{\\tilde{y}}\\ \\ (\\tilde{y}-\\hat{y})^{T}\\mathbf{\\Sigma}(\\tilde{y}-\\hat{y})\\qquad\\textrm{subject to}\\ \\ \\tilde{y}=\\mathbf{S}\\mathbf{G}\\tilde{y}," } }, { "type": "Paragraph", "content": [ "where ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚺∈ℝn×n", "meta": { "altText": "\\mathbf{\\Sigma}\\in\\mathbb{R}^{n\\times n}" } }, " is the covariance matrix of the coherency errors\n", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y~−y^", "meta": { "altText": "\\tilde{y}-\\hat{y}" } }, ", and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝐆=[𝟎n×(n−m)𝐈m]", "meta": { "altText": "\\mathbf{G}=[\\mathbf{0}_{n\\times(n-m)}\\quad\\mathbf{I}_{m}]" } }, ".\nIf ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚺", "meta": { "altText": "\\mathbf{\\Sigma}" } }, " is assumed to be known, the problem has a closed form solution, given by" ] }, { "type": "MathBlock", "id": "S1.E2", "label": "(2)", "mathLanguage": "mathml", "text": "y~=𝐒⁢(𝐒T⁢𝚺−1⁢𝐒)−1⁢𝐒T⁢𝚺−1⁢y^.", "meta": { "altText": "\\tilde{y}=\\mathbf{S}(\\mathbf{S}^{T}\\mathbf{\\Sigma}^{-1}\\mathbf{S})^{-1}\\mathbf{S}^{T}\\mathbf{\\Sigma}^{-1}\\hat{y}." } }, { "type": "Paragraph", "content": [ "In [", { "type": "Cite", "target": "bib-bib17", "content": [ "17" ] }, "] the authors showed that ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚺", "meta": { "altText": "\\mathbf{\\Sigma}" } }, " is not identifiable; instead, they propose\nthe trace minimization approach (MinT), which uses the covariance matrix of the base forecast\nerrors, i.e.," ] }, { "type": "MathBlock", "id": "S1.Ex3", "mathLanguage": "mathml", "text": "𝚺=V⁢[ε],", "meta": { "altText": "\\mathbf{\\Sigma}=V[\\varepsilon]," } }, { "type": "Paragraph", "content": [ "where ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "ε=y−y^", "meta": { "altText": "\\varepsilon=y-\\hat{y}" } }, " with ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y∈ℝn", "meta": { "altText": "y\\in\\mathbb{R}^{n}" } }, " being the observations arranged similarly to the base forecasts\nin equation (", { "type": "Cite", "target": "S1-E1", "content": [ "1" ] }, ").\nThis approach has shown great promise in increasing forecast accuracy across all levels of the hierarchy. MinT-Shrink,\nalso proposed by [", { "type": "Cite", "target": "bib-bib17", "content": [ "17" ] }, "], and extended to temporal hierarchies by [", { "type": "Cite", "target": "bib-bib14", "content": [ "14" ] }, "], is an\nextension of the MinT approach which applies a shrinkage estimator to the correlation." ] }, { "type": "Paragraph", "content": [ "The method estimates the covariance ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚺shrink", "meta": { "altText": "\\mathbf{\\Sigma}^{\\mathrm{shrink}}" } }, " from the correlation matrix of the base\nforecast errors ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝐑", "meta": { "altText": "\\mathbf{R}" } }, " using the shrinkage parameter ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "λ", "meta": { "altText": "\\lambda" } }, ":" ] }, { "type": "MathBlock", "id": "S1.E3", "label": "(3)", "mathLanguage": "mathml", "text": "𝐑shrink=(1−λ)⁢𝐑+λ⁢𝐈n,", "meta": { "altText": "\\displaystyle\\mathbf{R}^{\\mathrm{shrink}}=(1-\\lambda)\\mathbf{R}+\\lambda\\mathbf{I}_{n}," } }, { "type": "MathBlock", "id": "S1.E4", "label": "(4)", "mathLanguage": "mathml", "text": "𝚺shrink=𝚲12⁢𝐑shrink⁢𝚲12,", "meta": { "altText": "\\displaystyle\\mathbf{\\Sigma}^{\\mathrm{shrink}}=\\mathbf{\\Lambda}^{\\frac{1}{2}}\\mathbf{R}^{\\mathrm{shrink}}\\mathbf{\\Lambda}^{\\frac{1}{2}}," } }, { "type": "Paragraph", "content": [ "where ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲", "meta": { "altText": "\\mathbf{\\Lambda}" } }, " is the so-called hierarchical variance, which is a diagonal matrix consisting of the\nempirical variances of the base forecast errors at each level. So, for the example hierarchy in Figure ", { "type": "Cite", "target": "S1-F1", "content": [ "1" ] }, ",\nthe hierarchical\nvariance is" ] }, { "type": "MathBlock", "id": "S1.Ex4", "mathLanguage": "mathml", "text": "𝚲=[(σ24⁢h24⁢h)200000(σ24⁢h12⁢h)200000(σ12⁢h12⁢h)200000(σ24⁢h6⁢h)200000(σ18⁢h6⁢h)200000(σ12⁢h6⁢h)200000(σ6⁢h6⁢h)2].", "meta": { "altText": "\\mathbf{\\Lambda}=\\begin{bmatrix}\\bigl(\\sigma_{24\\mathrm{h}}^{24\\mathrm{h}}\\bigr)^{2}\\!\\!\\!\\!&0&0&0&0&0&\\!\\!\\!\\!0\\\\\n0\\!\\!\\!\\!&\\!\\!\\!\\!\\bigl(\\sigma_{24\\mathrm{h}}^{12\\mathrm{h}}\\bigr)^{2}\\!\\!\\!\\!&0&0&0&0&\\!\\!\\!\\!0\\\\\n0\\!\\!\\!\\!&0&\\!\\!\\!\\!\\bigl(\\sigma_{12\\mathrm{h}}^{12\\mathrm{h}}\\bigr)^{2}\\!\\!\\!\\!&0&0&0&\\!\\!\\!\\!0\\\\\n0\\!\\!\\!\\!&0&0&\\!\\!\\!\\!\\bigl(\\sigma_{24\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\!\\!\\!\\!&0&0&\\!\\!\\!\\!0\\\\\n0\\!\\!\\!\\!&0&0&0&\\!\\!\\!\\!\\bigl(\\sigma_{18\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\!\\!\\!\\!&0&\\!\\!\\!\\!0\\\\\n0\\!\\!\\!\\!&0&0&0&0&\\!\\!\\!\\!\\bigl(\\sigma_{12\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\!\\!\\!\\!&\\!\\!\\!\\!0\\\\\n0\\!\\!\\!\\!&0&0&0&0&0&\\!\\!\\!\\!\\bigl(\\sigma_{6\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\end{bmatrix}." } }, { "type": "Paragraph", "content": [ "It is usually assumed that both ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲", "meta": { "altText": "\\mathbf{\\Lambda}" } }, " and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝐑", "meta": { "altText": "\\mathbf{R}" } }, " are either constant or slowly varying\n(see, e.g., [", { "type": "Cite", "target": "bib-bib3", "content": [ "3" ] }, "]) in time. However, as mentioned in the introduction, this might not be a\nreasonable assumption for wind power forecasts, since the uncertainty is directly coupled to the wind power generation\nthrough the power curve, as discussed by, e.g., [", { "type": "Cite", "target": "bib-bib10", "content": [ "10" ] }, "]. A similar relation between the uncertainty and the\nprediction is probably also needed in other cases, like for solar power forecasting." ] }, { "type": "Heading", "id": "S2", "depth": 1, "content": [ "2 Wind power data" ] }, { "type": "Paragraph", "id": "S2.p1", "content": [ "The hierarchy that will be explored in this study has eight layers ranging from the 1-hour resolution up to a 24-hour\nresolution, i.e., the layers will have resolutions 1, 2, 3, 4, 6, 8, 12, and 24 hours. To perform the reconciliation,\nbase forecasts need to be obtained for each temporal resolution in the hierarchy." ] }, { "type": "Paragraph", "id": "S2.p2", "content": [ "The data that these forecasts should be generated from consists of hourly measurements of wind power production from\non-shore wind turbines, as well as numerical weather predictions (NWPs) of wind speed and wind direction. These data\nare available for 15 areas of the DK1 region in western Denmark (Jutland and Funen), as seen in Figure ", { "type": "Cite", "target": "S2-F2", "content": [ "2" ] }, ".\nThe areas have very different capacities for power production, ranging from as low as ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "39.3", "meta": { "altText": "39.3" } }, " MW to\n", { "type": "MathFragment", "mathLanguage": "mathml", "text": "746.7", "meta": { "altText": "746.7" } }, " MW, and each of the areas will be forecasted and reconciled individually using temporal hierarchies. The same\nregion has also been examined by [", { "type": "Cite", "target": "bib-bib6", "content": [ "6" ] }, "] using a spatial hierarchy." ] }, { "type": "Paragraph", "id": "S2.p3", "content": [ "Data were available from the beginning of January 2018 and through December 2019, and hence the data were split such\nthat 2018 data were used for training base forecast models and estimating shrinkage parameters, while 2019 data are\nused for testing the reconciliation out of sample.\nOut of sample testing will be based on a rolling window containing 12 months of previous measurements and forecasts.\nThis window rolls forward in time, one day at a time, as the data would become available. Each day the coefficients of\nthe base forecast models are re-estimated, and new base forecasts for the following day are produced and reconciled.\nThis mimics how operational wind power forecasting takes place." ] }, { "type": "Figure", "id": "S2-F2", "caption": [ { "type": "Paragraph", "content": [ "Map of the 15 subareas of the DK1 price area in western Denmark. Coastline in blue." ] } ], "licenses": [ { "type": "CreativeWork", "url": "https://creativecommons.org/licenses/by/4.0/legalcode", "content": [ { "type": "Paragraph", "content": [] } ] } ], "content": [ { "type": "ImageObject", "contentUrl": "dk1map_v3n.svg", "mediaType": "image/svg+xml", "meta": { "inline": false } } ] }, { "type": "Table", "id": "S2-T1", "caption": [ { "type": "Paragraph", "content": [ "Rated power production for each of the 15 subareas in the DK1 price area given in MW." ] } ], "rows": [ { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "Area" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "Rated power [MW]" ] } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "1" ] }, { "type": "TableCell", "content": [ "428.1" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "2" ] }, { "type": "TableCell", "content": [ "400.9" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "3" ] }, { "type": "TableCell", "content": [ "253.3" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "4" ] }, { "type": "TableCell", "content": [ "360.8" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "5" ] }, { "type": "TableCell", "content": [ "746.7" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "6" ] }, { "type": "TableCell", "content": [ "328.0" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "7" ] }, { "type": "TableCell", "content": [ "86.9" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "8" ] }, { "type": "TableCell", "content": [ "100.4" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "9" ] }, { "type": "TableCell", "content": [ "39.3" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "10" ] }, { "type": "TableCell", "content": [ "201.6" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "11" ] }, { "type": "TableCell", "content": [ "219.7" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "12" ] }, { "type": "TableCell", "content": [ "97.3" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "13" ] }, { "type": "TableCell", "content": [ "194.4" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "14" ] }, { "type": "TableCell", "content": [ "230.2" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "15" ] }, { "type": "TableCell", "content": [ "75.6" ] } ] } ] }, { "type": "Paragraph", "id": "S2.p4", "content": [ "From each of the areas depicted, hourly measurements of wind power production, as well as forecasts of wind speed and\ndirection, are available. The power production data is courtesy of Energinet, the TSO that operates the Danish\nhigh-voltage transmission grids. The weather forecasts were made by the Danish Meteorological Institute (DMI)." ] }, { "type": "Paragraph", "id": "S2.p5", "content": [ "Commercial state-of-the-art forecasts of hourly wind power production were available for the 1-hour resolution\n(courtesy of ENFOR A/S). These were produced using their wind power forecasting tool\nWindFor,", { "type": "Note", "id": "idm1048", "noteType": "Footnote", "content": [ { "type": "Paragraph", "id": "footnote2", "content": [ { "type": "Link", "target": "https://enfor.dk/services/windfor/", "content": [ "https://enfor.dk/services/windfor/" ] } ] } ] }, " and will constitute the lowest level of the temporal hierarchy." ] }, { "type": "Paragraph", "id": "S2.p6", "content": [ "Table ", { "type": "Cite", "target": "S2-T2", "content": [ "2" ] }, " below summarizes which data are available along with the notations in this study." ] }, { "type": "Table", "id": "S2-T2", "caption": [ { "type": "Paragraph", "content": [ "Data descriptions for variables and the used units. These are all on 1-hour time resolution for each of the 15 areas." ] } ], "rows": [ { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "Variable" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "Unit" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "Description" ] } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "y", "meta": { "altText": "y" } } ] }, { "type": "TableCell", "content": [ "kWh" ] }, { "type": "TableCell", "content": [ "Measurements of power production" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "y^", "meta": { "altText": "\\hat{y}" } } ] }, { "type": "TableCell", "content": [ "kWh" ] }, { "type": "TableCell", "content": [ "Forecasts of power production" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "WSpd^", "meta": { "altText": "\\widehat{\\mathrm{WSpd}}" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "m/s", "meta": { "altText": "\\mathrm{m}/\\mathrm{s}" } } ] }, { "type": "TableCell", "content": [ "Forecasted wind speed" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "WDir^", "meta": { "altText": "\\widehat{\\mathrm{WDir}}" } } ] }, { "type": "TableCell", "content": [ "degrees" ] }, { "type": "TableCell", "content": [ "Forecasted wind direction" ] } ] } ] }, { "type": "Paragraph", "id": "S2.p7", "content": [ "The measurements ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y", "meta": { "altText": "y" } }, " are settlement measurements, meaning that they are generally available 8–10 days after the\nmeasured production has occurred. This fact will have to be accounted for when modelling the power production, since\nnot having real-time data means that some models, such as the classic ARIMA models, are not applicable. The\nmeasurements also include periods of missing data. For shorter periods (", { "type": "MathFragment", "mathLanguage": "mathml", "text": "≤6", "meta": { "altText": "\\leq 6" } }, " hours), measurements were linearly\ninterpolated, while longer periods were left missing, meaning that models could not be updated here." ] }, { "type": "Paragraph", "id": "S2.p8", "content": [ "The forecasts in this study are generated every night at midnight (00:00) and cover the following 24-hour period. The\nforecasts of wind speed and direction are generated every six hours and are available a few hours later. To generate\nforecasts at midnight, the NWPs for the following day were taken from the weather forecast generated at 18:00. This\nalso means that the forecast horizon follows the time of day." ] }, { "type": "Figure", "id": "S2-F3", "caption": [ { "type": "Paragraph", "content": [ "Left: Time series plot of the wind power production in area no. 3 in northern Jutland from\nJanuary 2018. Right: Relation between forecasted wind speed for all horizons and their corresponding normalized power\nproductions in area no. 3." ] } ], "content": [ { "type": "ImageObject", "contentUrl": "172-f3.png", "mediaType": "image/png", "meta": { "inline": false } } ] }, { "type": "Paragraph", "id": "S2.p9", "content": [ "The time series in Figure ", { "type": "Cite", "target": "S2-F3", "content": [ "3" ] }, " shows some of the volatility in the data, as changes from low\nproduction to high production and vice versa can happen within a couple of hours." ] }, { "type": "Paragraph", "id": "S2.p10", "content": [ "The relation between wind speed and power production seen in Figure ", { "type": "Cite", "target": "S2-F3", "content": [ "3" ] }, " shows the classic power\ncurve structure, with power production increasing slowly at first until about 4 m/s, then rapidly increasing above the\n4 m/s mark, and flattening out when the forecasted wind speed approaches 8 m/s. This power curve structure is also what\nshould be mimicked when modelling the individual layers of the hierarchy." ] }, { "type": "Paragraph", "id": "S2.p11", "content": [ "Figure ", { "type": "Cite", "target": "S2-F3", "content": [ "3" ] }, " also illustrates how the variance in power production depends on the wind speed, and\nhence on the predicted wind power generation. This means that the wind power production cannot be assumed to be\nGaussian, since the variance would then be independent of the power production." ] }, { "type": "Paragraph", "id": "S2.p12", "content": [ "Lastly, it is worth noting that the power production does not reach the rated production in the available data. This is\nlikely due to the fact that the turbines within the areas are positioned differently in the terrain, making it very\nrare for all of them to generate according to their rated power simultaneously." ] }, { "type": "Heading", "id": "S2.SS1", "depth": 2, "content": [ "2.1 Aggregation" ] }, { "type": "Paragraph", "id": "S2.SS1.p1", "content": [ "In order to construct forecasts for the other levels of the hierarchy, the data had to be aggregated up to the 24-hour\nresolution. Thus, datasets of power production and weather predictions need to be constructed for the 2, 3, 4, 6, 8,\n12, and 24-hour resolutions.\nThis is done by aggregating the 1-hour measurements and NWPs. Power production measurements are aggregated by simply\nadding up measurements, while wind speed and direction forecasts are aggregated by averaging the forecasts. For wind\nspeed, this is a simple arithmetic average, while for the direction, the average has to be computed using the circular\nmean. This implies converting the angular directions to points on the unit circle, taking the arithmetic mean of the\npoints, with the resulting angle being the circular mean." ] }, { "type": "Paragraph", "id": "S2.SS1.p2", "content": [ "At the 1-hour level, the 12 months of training data contain ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "8760", "meta": { "altText": "8760" } }, " measurements of power production and as many\nforecasted values of wind speed and direction. As a result of the aggregation, the 2-hour level has ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "4380", "meta": { "altText": "4380" } }, " data points,\nthe 3-hour level has ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "2920", "meta": { "altText": "2920" } }, " data points, and so on, up to the 24-hour level, which has ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "365", "meta": { "altText": "365" } }, " data points. This\nreduction in data quantity puts a limit on how well the higher levels can be modelled. Consequently, if a larger\ntemporal hierarchy should be attempted, more data would be required." ] }, { "type": "Heading", "id": "S3", "depth": 1, "content": [ "3 Base forecasts" ] }, { "type": "Paragraph", "id": "S3.p1", "content": [ "When modelling the levels above the 1-hour resolution, i.e., 2, 3, 4, 6, 8, 12, and 24 hours, it was decided to use\nbeta regression models. These purely statistical models are chosen mainly to address the challenge of\npower measurements not being available in real time.\nHowever, this is also a well-established method of predicting wind energy, see, e.g., [", { "type": "Cite", "target": "bib-bib15", "content": [ "15" ] }, "]. Additionally,\nchoosing such simple models can also show if improvements to already high-quality commercial forecasts are possible\nusing forecast reconciliation, even without employing state-of-the-art forecasts for the other levels.\n" ] }, { "type": "Paragraph", "id": "S3.p2", "content": [ "The beta regression models used in this study are constructed using the betareg package in R [", { "type": "Cite", "target": "bib-bib4", "content": [ "4" ] }, "].\nThe following brief introduction leans heavily on the documentation for the package, so for a more in-depth\nwalk-through of these methods, see the documentation by [", { "type": "Cite", "target": "bib-bib4", "content": [ "4" ] }, "]." ] }, { "type": "Paragraph", "content": [ "We model the normalized power production as a beta distribution ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "Yt∼ℬ⁢(μt,ϕt)", "meta": { "altText": "Y_{t}\\sim\\mathcal{B}(\\mu_{t},\\phi_{t})" } }, ", where\n", { "type": "MathFragment", "mathLanguage": "mathml", "text": "μt∈(0,1)", "meta": { "altText": "\\mu_{t}\\in(0,1)" } }, " is the mean value (", { "type": "MathFragment", "mathLanguage": "mathml", "text": "E⁢[Yt]=μt", "meta": { "altText": "E[Y_{t}]=\\mu_{t}" } }, ") and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "0\" display=\"inline\">ϕt>0", "meta": { "altText": "\\phi_{t}>0" } }, " is the precision. Therefore, the variance of the power\nproduction at time ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, " can be expressed as" ] }, { "type": "MathBlock", "id": "S3.E5", "label": "(5)", "mathLanguage": "mathml", "text": "V⁢[Yt]=μt⁢(1−μt)1+ϕt.", "meta": { "altText": "V[Y_{t}]=\\frac{\\mu_{t}(1-\\mu_{t})}{1+\\phi_{t}}." } }, { "type": "Paragraph", "content": [ "In order to be able to introduce explanatory variables the parameters ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "μt", "meta": { "altText": "\\mu_{t}" } }, " and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "ϕt", "meta": { "altText": "\\phi_{t}" } }, " are modelled\nindividually by" ] }, { "type": "MathBlock", "id": "S3.Ex5", "mathLanguage": "mathml", "text": "g1⁢(μt)=xtT⁢β,g2⁢(ϕt)=ztT⁢γ.", "meta": { "altText": "g_{1}(\\mu_{t})=x_{t}^{T}\\beta,g_{2}(\\phi_{t})=z_{t}^{T}\\gamma." } }, { "type": "Paragraph", "content": [ "Here ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "β", "meta": { "altText": "\\beta" } }, " and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "γ", "meta": { "altText": "\\gamma" } }, " are regression coefficients, and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "xt", "meta": { "altText": "x_{t}" } }, " and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "zt", "meta": { "altText": "z_{t}" } }, " are regressors; ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "g1", "meta": { "altText": "g_{1}" } }, " and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "g2", "meta": { "altText": "g_{2}" } }, " are link\nfunctions [", { "type": "Cite", "target": "bib-bib11", "content": [ "11" ] }, "], where ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "g1", "meta": { "altText": "g_{1}" } }, " was chosen to be the logit function and\n", { "type": "MathFragment", "mathLanguage": "mathml", "text": "g2", "meta": { "altText": "g_{2}" } }, " to be the natural logarithm.\nGiven the assumed density, the coefficients ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "β", "meta": { "altText": "\\beta" } }, " and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "γ", "meta": { "altText": "\\gamma" } }, " are estimated by maximum likelihood estimation." ] }, { "type": "Paragraph", "content": [ "When constructing the models for the hierarchy, a method of forward selection was used based on the Wald test with\nsignificance level ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "p=0.05", "meta": { "altText": "p=0.05" } }, ".\nFor modelling the mean ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "μt", "meta": { "altText": "\\mu_{t}" } }, ", the initial model consisted only of the forecasted wind speed. The two following\nmodelling steps then added the second and third power of the wind speed. The fourth power did not provide a significant\nimprovement and was thus dropped. Next, the contribution of the wind direction was considered by adding the sine and\ncosine of the wind direction. However, only the cosine was significant, resulting in the following final model for the\nconditional mean:" ] }, { "type": "MathBlock", "id": "S3.Ex6", "mathLanguage": "mathml", "text": "g1⁢(μt+h|t)=β0+β1⁢WSpd^t+h|t+β2⁢WSpd^t+h|t2", "meta": { "altText": "\\displaystyle g_{1}(\\mu_{t+h|t})=\\beta_{0}+\\beta_{1}\\widehat{\\mathrm{WSpd}}_{t+h|t}+\\beta_{2}\\widehat{\\mathrm{WSpd}}_{t+h|t}^{2}" } }, { "type": "MathBlock", "id": "S3.Ex7", "mathLanguage": "mathml", "text": "+β3⁢WSpd^t+h|t3+β4⁢cos⁡(WDir^t+h|t),", "meta": { "altText": "\\displaystyle+\\beta_{3}\\widehat{\\mathrm{WSpd}}_{t+h|t}^{3}+\\beta_{4}\\cos\\bigl(\\widehat{\\mathrm{WDir}}_{t+h|t}\\bigr)," } }, { "type": "Paragraph", "content": [ "where\n", { "type": "MathFragment", "mathLanguage": "mathml", "text": "WSpd^t+h|t", "meta": { "altText": "\\widehat{\\mathrm{WSpd}}_{t+h|t}" } }, " is the forecasted wind speed from the NWP available at time ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, " for horizon ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "h", "meta": { "altText": "h" } }, "\nand ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "WDir^t+h|t", "meta": { "altText": "\\widehat{\\mathrm{WDir}}_{t+h|t}" } }, " is similarly the forecasted wind direction. Since the forecasts are produced daily\nat midnight, ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, " changes in steps of 24 hours, i.e., ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t=0,24,48,72,…", "meta": { "altText": "t=0,24,48,72,\\ldots" } }, "." ] }, { "type": "Paragraph", "content": [ "When modelling the precision ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "ϕ", "meta": { "altText": "\\phi" } }, ", a very similar procedure was used. The wind speed terms were significant up to the\nsecond order, and only the sine of the wind direction was significant. Additionally, the prediction horizon of the\nweather forecast was added as a regressor. This was added to address the general tendency that the uncertainty is\nincreasing with the prediction horizon. In the end, the precision is modelled as" ] }, { "type": "MathBlock", "id": "S3.E6", "label": "(6)", "mathLanguage": "mathml", "text": "g2⁢(ϕt+h|t)=γ0+γ1⁢WSpd^t+h|t+γ2⁢WSpd^t+h|t2", "meta": { "altText": "g_{2}(\\phi_{t+h|t})=\\gamma_{0}+\\gamma_{1}\\widehat{\\mathrm{WSpd}}_{t+h|t}+\\gamma_{2}\\widehat{\\mathrm{WSpd}}_{t+h|t}^{2}" } }, { "type": "MathBlock", "id": "S3.Ex8", "mathLanguage": "mathml", "text": "+γ3⁢sin⁡(WDir^t+h|t)+γ4⁢h,", "meta": { "altText": "{}+\\gamma_{3}\\sin\\bigl(\\widehat{\\mathrm{WDir}}_{t+h|t}\\bigr)+\\gamma_{4}h," } }, { "type": "Paragraph", "content": [ "where ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "h", "meta": { "altText": "h" } }, " is the prediction horizon of the weather forecast." ] }, { "type": "Figure", "id": "S3-F4", "caption": [ { "type": "Paragraph", "content": [ "Top: Power curve resulting from fitting the beta regression model to the aggregated 2-hour data for area\nno. 3. Bottom: Standard deviation of fitted values from beta regression model on aggregated 2-hour data for area no. 3,\ngiven the wind direction is fixed in a northerly direction." ] } ], "content": [ { "type": "ImageObject", "contentUrl": "172-f4.png", "mediaType": "image/png", "meta": { "inline": false } } ] }, { "type": "Paragraph", "id": "S3.p7", "content": [ "The power curve (Figure ", { "type": "Cite", "target": "S3-F4", "content": [ "4" ] }, ") resulting from fitting the model to the 2018 training set seems\nto fit well with what was expected from the data (Figure ", { "type": "Cite", "target": "S2-F3", "content": [ "3" ] }, "). However, it seems that the wind\ndirection plays a very small role, as the four wind directions (N, E, S, W) are indistinguishable\n(Figure ", { "type": "Cite", "target": "S3-F4", "content": [ "4" ] }, ", top)." ] }, { "type": "Paragraph", "id": "S3.p8", "content": [ "With forecasted weather variables, the uncertainty of the NWP is translated to power by the non-linear power curve. As\nnicely illustrated by [", { "type": "Cite", "target": "bib-bib10", "content": [ "10" ] }, "], this explains the form of the standard deviation in Figure ", { "type": "Cite", "target": "S3-F4", "content": [ "4" ] }, "\n(bottom), and this is also in line with what was seen by [", { "type": "Cite", "target": "bib-bib12", "content": [ "12" ] }, "]. Furthermore, the standard deviation\nclearly depends on the prediction horizon." ] }, { "type": "Paragraph", "content": [ "Examining the residuals of the fitted models for the area on each resolution, the accuracy can be assessed by the\n", { "type": "Emphasis", "content": [ "root-mean-square error" ] }, " (RMSE), which is\ncomputed as" ] }, { "type": "MathBlock", "id": "S3.Ex9", "mathLanguage": "mathml", "text": "RMSEarea,l=∑t=1Nl(yl,t−y^l,t)2Nl,", "meta": { "altText": "\\mathrm{RMSE}_{\\mathrm{area},l}=\\sqrt{\\sum_{t=1}^{N_{l}}\\frac{(y_{l,t}-\\hat{y}_{l,t})^{2}}{N_{l}}}," } }, { "type": "Paragraph", "content": [ "where ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "yt", "meta": { "altText": "y_{t}" } }, " are the wind power measurements at time ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, ", ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y^t", "meta": { "altText": "\\hat{y}_{t}" } }, " are the predicted values at time ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, ", and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "Nl", "meta": { "altText": "N_{l}" } }, " is\nthe amount of data points at resolution ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "l", "meta": { "altText": "l" } }, ", i.e. for the 2018 training set ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "N1⁢h=8760", "meta": { "altText": "N_{1\\mathrm{h}}=\\nobreak 8760" } }, ", ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "N2⁢h=4380", "meta": { "altText": "N_{2\\mathrm{h}}=\\nobreak 4380" } }, " etc." ] }, { "type": "Table", "id": "S3-T3", "caption": [ { "type": "Paragraph", "content": [ "Accuracy of the different models for area no. 3 on the 2018 training data. 1-hour predictions were made by\na state-of-the-art commercial model, while those for 2 hours and beyond are from the beta regressions." ] } ], "rows": [ { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "Model resolution" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "RMSE[MW]" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "% of rated power" ] } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "1 hour" ] }, { "type": "TableCell", "content": [ "15.50" ] }, { "type": "TableCell", "content": [ "6.12" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "2 hours" ] }, { "type": "TableCell", "content": [ "33.39" ] }, { "type": "TableCell", "content": [ "6.59" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "3 hours" ] }, { "type": "TableCell", "content": [ "48.04" ] }, { "type": "TableCell", "content": [ "6.32" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "4 hours" ] }, { "type": "TableCell", "content": [ "62.26" ] }, { "type": "TableCell", "content": [ "6.14" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "6 hours" ] }, { "type": "TableCell", "content": [ "88.89" ] }, { "type": "TableCell", "content": [ "5.85" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "8 hours" ] }, { "type": "TableCell", "content": [ "112.58" ] }, { "type": "TableCell", "content": [ "5.55" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "12 hours" ] }, { "type": "TableCell", "content": [ "155.79" ] }, { "type": "TableCell", "content": [ "5.12" ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "24 hours" ] }, { "type": "TableCell", "content": [ "265.12" ] }, { "type": "TableCell", "content": [ "4.36" ] } ] } ] }, { "type": "Paragraph", "id": "S3.p10", "content": [ "While the RMSE values as seen in Table ", { "type": "Cite", "target": "S3-T3", "content": [ "3" ] }, " increase gradually with the temporal resolution, it is\nseen that relative to the rated power all forecasts perform well. It should, however, be noted that this is not a\nfair comparison between the 1-hour forecast and the other resolutions, as the forecasts on the 1-hour resolution were\nmade out of sample in real time, while the forecasts on the other resolutions are in sample." ] }, { "type": "Heading", "id": "S4", "depth": 1, "content": [ "4 Forecast reconciliation methodology" ] }, { "type": "Paragraph", "id": "S4.p1", "content": [ "With the base forecasts that have been generated, the forecasts in the temporal hierarchy can be reconciled. With the\nMinT-Shrink method, this implies computing the hierarchical variance ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲", "meta": { "altText": "\\mathbf{\\Lambda}" } }, " and the correlation\nmatrix ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝐑", "meta": { "altText": "\\mathbf{R}" } }, ", inserting them in the equations (", { "type": "Cite", "target": "S1-E3", "content": [ "3" ] }, ") and (", { "type": "Cite", "target": "S1-E4", "content": [ "4" ] }, "), and finally\ncomputing the reconciled forecast by equation (", { "type": "Cite", "target": "S1-E2", "content": [ "2" ] }, ")." ] }, { "type": "Paragraph", "id": "S4.p2", "content": [ "This on its own is likely to improve the accuracy of the base forecasts, as seen in, e.g., [", { "type": "Cite", "target": "bib-bib5", "content": [ "5" ] }, "] for wind\npower. However, we propose an alternative which better addresses some of the inherent challenges calling for a\nfunctional relationship between the uncertainty (conditional variance) and the predicted wind power generation\n(conditional mean)." ] }, { "type": "Paragraph", "id": "S4.p3", "content": [ "One such challenge is that the boundedness of wind power production should be accounted for, since wind power is\ninherently bounded both from below at zero and from above at the rated power production. This means that the covariance\nwill depend on the forecasted power production. This dependence is not accounted for in the MinT-Shrink method, where\nthe covariance is assumed to be constant." ] }, { "type": "Paragraph", "id": "S4.p4", "content": [ "To address this challenge, we propose a method for linking the hierarchical variance ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲", "meta": { "altText": "\\mathbf{\\Lambda}" } }, " to the\npredicted power production." ] }, { "type": "Paragraph", "content": [ "In each time step ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, ", ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲t", "meta": { "altText": "\\mathbf{\\Lambda}_{t}" } }, " will be set based on a parametric estimate of the prediction uncertainty.\nIn the case of the example hierarchy from Figure ", { "type": "Cite", "target": "S1-F1", "content": [ "1" ] }, ", at time ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, ",\n", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲t", "meta": { "altText": "\\mathbf{\\Lambda}_{t}" } }, " becomes" ] }, { "type": "MathBlock", "id": "S4.Ex10", "mathLanguage": "mathml", "text": "𝚲t=diag⁢([(σ^t+24⁢h24⁢h)2(σ^t+24⁢h12⁢h)2(σ^t+12⁢h12⁢h)2(σ^t+24⁢h6⁢h)2(σ^t+18⁢h6⁢h)2(σ^t+12⁢h6⁢h)2(σ^t+6⁢h6⁢h)2]).", "meta": { "altText": "\\mathbf{\\Lambda}_{t}={\\mathrm{diag}}\\left(\\begin{bmatrix}\\bigl(\\hat{\\sigma}_{t+24\\mathrm{h}}^{24\\mathrm{h}}\\bigr)^{2}\\\\[2.84526pt]\n\\bigl(\\hat{\\sigma}_{t+24\\mathrm{h}}^{12\\mathrm{h}}\\bigr)^{2}\\\\[2.84526pt]\n\\bigl(\\hat{\\sigma}_{t+12\\mathrm{h}}^{12\\mathrm{h}}\\bigr)^{2}\\\\[2.84526pt]\n\\bigl(\\hat{\\sigma}_{t+24\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\\\[2.84526pt]\n\\bigl(\\hat{\\sigma}_{t+18\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\\\[2.84526pt]\n\\bigl(\\hat{\\sigma}_{t+12\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\\\[2.84526pt]\n\\bigl(\\hat{\\sigma}_{t+6\\mathrm{h}}^{6\\mathrm{h}}\\bigr)^{2}\\end{bmatrix}\\right)." } }, { "type": "Paragraph", "content": [ "This forecast-dependent hierarchical variance ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲t", "meta": { "altText": "\\mathbf{\\Lambda}_{t}" } }, " then replaces the constant hierarchical\nvariance ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝚲", "meta": { "altText": "\\mathbf{\\Lambda}" } }, ". The shrinkage equations thus\nbecome" ] }, { "type": "MathBlock", "id": "S4.Ex11", "mathLanguage": "mathml", "text": "𝐑Var=(1−λVar)⁢𝐑+λVar⁢𝐈n,𝚺tVar=𝚲t12⁢𝐑Var⁢𝚲t12,", "meta": { "altText": "\\mathbf{R}^{\\mathrm{Var}}=(1-\\lambda^{\\mathrm{Var}})\\mathbf{R}+\\lambda^{\\mathrm{Var}}\\mathbf{I}_{n},\\qquad\\mathbf{\\Sigma}^{\\mathrm{Var}}_{t}=\\mathbf{\\Lambda}_{t}^{\\frac{1}{2}}\\mathbf{R}^{\\mathrm{Var}}\\mathbf{\\Lambda}_{t}^{\\frac{1}{2}}," } }, { "type": "Paragraph", "content": [ "which\nwill be denoted MinT-Var." ] }, { "type": "Paragraph", "id": "S4.p6", "content": [ "This approach adapts the variance to the forecast, but not the correlation, which is still estimated from the response\nresiduals. The next step is to allow the correlation structure to change with the forecast. One way of doing this\nis to compute the correlation from Pearson residuals instead, which are normalized by the forecast-dependent variance." ] }, { "type": "Paragraph", "id": "S4.p7", "content": [ "Additionally, it is well-known from the theory of generalized linear models that Pearson residuals will be closer to\nnormal distributions than the original response residuals, see e.g. [", { "type": "Cite", "target": "bib-bib11", "content": [ "11" ] }, "]." ] }, { "type": "Paragraph", "content": [ "For a set of forecasts at time ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, ", the Pearson residuals may be\nexpressed as" ] }, { "type": "MathBlock", "id": "S4.Ex12", "mathLanguage": "mathml", "text": "εtp=𝚲t−1/2⁢(yt−yt^)∼N⁢(0,𝚺p).", "meta": { "altText": "\\varepsilon_{t}^{p}=\\mathbf{\\Lambda}_{t}^{-1/2}(y_{t}-\\hat{y_{t}})\\sim N(0,\\mathbf{\\Sigma}^{p})." } }, { "type": "Paragraph", "content": [ "The second proposal is thus to use both the correlation of the Pearson residuals ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "Rp", "meta": { "altText": "R^{p}" } }, " and a forecast-dependent\nvariance, i.e." ] }, { "type": "MathBlock", "id": "S4.Ex13", "mathLanguage": "mathml", "text": "𝐑PVar=(1−λPVar)⁢𝐑p+λPVar⁢𝐈n,𝚺tPVar=𝚲t12⁢𝐑PVar⁢𝚲t12,", "meta": { "altText": "\\mathbf{R}^{\\mathrm{PVar}}=(1-\\lambda^{\\mathrm{PVar}})\\mathbf{R}^{p}+\\lambda^{\\mathrm{PVar}}\\mathbf{I}_{n},\\qquad\\mathbf{\\Sigma}^{\\mathrm{PVar}}_{t}=\\mathbf{\\Lambda}_{t}^{\\frac{1}{2}}\\mathbf{R}^{\\mathrm{PVar}}\\mathbf{\\Lambda}_{t}^{\\frac{1}{2}}," } }, { "type": "Paragraph", "content": [ "which will be denoted MinT-PVar." ] }, { "type": "Paragraph", "id": "S4.p10", "content": [ "These approaches lead to a flexible and more correct description of the uncertainty than the original method, while\nalso addressing the challenges imposed when reconciling wind power forecasts. The performance of the methods will,\nhowever, depend on how well the prediction uncertainty is estimated. This then raises the question of how one should\nestimate the prediction uncertainty." ] }, { "type": "Paragraph", "id": "S4.p11", "content": [ "In our case, where the beta regressions are used to generate base forecasts, the variance of the predictions is given\nby equation (", { "type": "Cite", "target": "S3-E5", "content": [ "5" ] }, "). Conversely, the variance on the 1-hour resolution is not known, since the\nunderlying model is not available. Therefore, for the 1-hour resolution the prediction uncertainty will have to be estimated\nby another method. Here, a well-known and effective way of doing this is through maximum likelihood estimation, which\nalso allows us to describe how the conditional variance depends on the conditional mean." ] }, { "type": "Paragraph", "content": [ "Assuming that residuals are normally distributed, a likelihood function can be established for each prediction\nhorizon ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "h", "meta": { "altText": "h" } }, ". For the set of ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "k", "meta": { "altText": "k" } }, " observations ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝒴t+h={yt−24⁢k+h,…,yt−24+h}", "meta": { "altText": "\\mathcal{Y}_{t+h}=\\{y_{t-24k+h},\\ldots,y_{t-24+h}\\}" } }, " and the\ncorresponding set of forecasts ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "𝒴^t+h={y^t−24⁢k+h,…,y^t−24+h}", "meta": { "altText": "\\hat{\\mathcal{Y}}_{t+h}=\\{\\hat{y}_{t-24k+h},\\ldots,\\hat{y}_{t-24+h}\\}" } }, ", the\nlog-likelihood will be\ncomputed as" ] }, { "type": "MathBlock", "id": "S4.Ex14", "mathLanguage": "mathml", "text": "ℓ⁢(σ^t+h2;𝒴^t+h,𝒴t+h)=∑i=1k−k2⁢log⁡(σ^t+h2)", "meta": { "altText": "\\displaystyle\\ell(\\hat{\\sigma}_{t+h}^{2};\\hat{\\mathcal{Y}}_{t+h},\\mathcal{Y}_{t+h})=\\sum_{i=1}^{k}-\\frac{k}{2}\\log(\\hat{\\sigma}_{t+h}^{2})" } }, { "type": "MathBlock", "id": "S4.Ex15", "mathLanguage": "mathml", "text": "+(−12⁢σ^t+h2⁢(yt−24⁢i+h−y^t−24⁢i+h)2).", "meta": { "altText": "\\displaystyle+\\Bigl(-\\frac{1}{2\\hat{\\sigma}_{t+h}^{2}}(y_{t-24i+h}-\\hat{y}_{t-24i+h})^{2}\\Bigr)." } }, { "type": "Paragraph", "content": [ "This requires a parametrization of the variances ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "σt2", "meta": { "altText": "\\sigma_{t}^{2}" } }, ", for which inspiration can be drawn from the beta\nregression models for the precision as well as the relation between variance and precision given in\nequation (", { "type": "Cite", "target": "S3-E5", "content": [ "5" ] }, ").\nThis results in the following parametrization of the\nvariance:" ] }, { "type": "MathBlock", "mathLanguage": "mathml", "text": "σ^t+h2=exp⁡(α0,h)+y^t+h⁢(1−y^t+h)1+exp⁡(α1,h+α2,h⁢WSpd^t+h+α3,h⁢sin⁡(WDir^t+h)).", "meta": { "altText": "\\hat{\\sigma}_{t+h}^{2}=\\exp(\\alpha_{0,h})\\\\\n{}+\\frac{\\hat{y}_{t+h}(1-\\hat{y}_{t+h})}{1+\\exp(\\alpha_{1,h}+\\alpha_{2,h}\\widehat{\\mathrm{WSpd}}_{t+h}+\\alpha_{3,h}\\sin(\\widehat{\\mathrm{WDir}}_{t+h}))}." } }, { "type": "Paragraph", "content": [ "This parametrization represents a slight simplification of the precision model (equation (", { "type": "Cite", "target": "S3-E6", "content": [ "6" ] }, ")),\nas the second-order term is omitted. This simplification is justified by examining the ", { "type": "Emphasis", "content": [ "Akaike information criterion" ] }, " (AIC)\nof the likelihood functions, with the full precision model and\nwith the simplified model both fitted on the training data. The simpler model had a lower AIC, and hence it was chosen\ninstead of the full precision model.\nFurthermore, this parametrization has an added constant term determined by ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "α0,h", "meta": { "altText": "\\alpha_{0,h}" } }, ". The addition of this term\nimproved the robustness of the optimization and lowered the AIC significantly. Moreover, it is in line with the\nbehaviour of the variance seen in Figure ", { "type": "Cite", "target": "S3-F4", "content": [ "4" ] }, " not going fully to zero when there is no wind.\nFigure ", { "type": "Cite", "target": "S4-F5", "content": [ "5" ] }, " shows contour plots of the estimated prediction uncertainty with parameters estimated in\nsample for area no. 3." ] }, { "type": "Figure", "id": "S4-F5", "caption": [ { "type": "Paragraph", "content": [ "Left: Contour plot of variance for wind speed and direction; 1-hour prediction horizon and normalized power production of 0.5. Right: Contour plot of variance for prediction horizon and power output; 8 m/s wind speed and ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "180∘", "meta": { "altText": "180^{\\circ}" } }, " wind direction." ] } ], "content": [ { "type": "ImageObject", "contentUrl": "172-f5.png", "mediaType": "image/png", "meta": { "inline": false } } ] }, { "type": "Paragraph", "id": "S4.p14", "content": [ "The uncertainty seems to be dependent linearly on the wind speed, while the uncertainty varies\nsinusoidally with the wind direction, topping at around ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "270∘", "meta": { "altText": "270^{\\circ}" } }, ", corresponding to a westerly wind direction. This is\nin line with the behaviour observed by, e.g., [", { "type": "Cite", "target": "bib-bib13", "content": [ "13" ] }, "]. The uncertainty related to power output quite\nclearly follows the polynomial relation from the variance model, with maximum at around ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "y^t+h=0.5", "meta": { "altText": "\\hat{y}_{t+h}=0.5" } }, ". For the\nprediction horizon (Figure ", { "type": "Cite", "target": "S4-F5", "content": [ "5" ] }, ", right column) the picture is less clear. The first few hours have quite\nlow uncertainty, increasing gradually before spiking around the 10-hour horizon, then dropping down again for a few\nhours before again spiking around the 19–21-hour mark." ] }, { "type": "Paragraph", "id": "S4.p15", "content": [ "To find the coefficients ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "αh", "meta": { "altText": "\\alpha_{h}" } }, " of the model, the log-likelihood is numerically optimized in each time step ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "t", "meta": { "altText": "t" } }, " for\neach horizon ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "h", "meta": { "altText": "h" } }, " using a window of past data. The resulting variance estimate can then be used in the reconciliation of\nthe forecast." ] }, { "type": "Paragraph", "id": "S4.p16", "content": [ "Lastly, to estimate the shrinkage intensities, i.e., ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "λ", "meta": { "altText": "\\lambda" } }, ", ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "λVar", "meta": { "altText": "\\lambda^{\\mathrm{Var}}" } }, ", ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "λPVar", "meta": { "altText": "\\lambda^{\\mathrm{PVar}}" } }, " we\nfollow [", { "type": "Cite", "target": "bib-bib14", "content": [ "14" ] }, "] and optimize them by doing a grid search, discretizing the regularization parameter on a\ngrid of ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "λ=0,0.01,0.02,…,1", "meta": { "altText": "\\lambda=0,0.01,0.02,\\ldots,1" } }, ". The parameters are then chosen based on the greatest in-sample improvements\nto RMSE for each area. The resulting optimal parameter is then used for the out-of-sample tests." ] }, { "type": "Heading", "id": "S5", "depth": 1, "content": [ "5 Results" ] }, { "type": "Paragraph", "id": "S5.p1", "content": [ "In this section, the out-of-sample performance of the three different methods for reconciling will be examined, those\nbeing the MinT-Shrink as proposed by [", { "type": "Cite", "target": "bib-bib17", "content": [ "17" ] }, "], and the two proposed in this study, MinT-Var and\nMinT-PVar.\nFor testing, a rolling window will move one day at a time. Each day, the coefficients of the base forecast models are\nre-estimated, and new base forecasts for the following day are produced. Then the variances of the base forecasts are\nestimated. In the case of the beta regressions, this is done directly using equation (", { "type": "Cite", "target": "S3-E5", "content": [ "5" ] }, "), and for\nthe commercial forecasts on the 1-hour resolution, maximum likelihood estimation is performed using the windowed data\nto estimate coefficients. Then the base forecasts are reconciled using each of the three methods. Lastly, the\ncovariances of the residuals are computed, so they can be used for the reconciliation of the following day’s data." ] }, { "type": "Paragraph", "content": [ "To evaluate the forecasts for individual areas, the RRMSE% (", { "type": "Emphasis", "content": [ "percentage relative root-mean-square error" ] }, ") is examined\nas a way to evaluate the performance and measure the accuracy of the proposed methods. For each area and for each\nresolution (", { "type": "MathFragment", "mathLanguage": "mathml", "text": "l", "meta": { "altText": "l" } }, "), the RRMSE% is\ngiven as" ] }, { "type": "MathBlock", "id": "S5.Ex16", "mathLanguage": "mathml", "text": "RRMSE%area,l=100%⁢(RRMSEarea,lRRMSEarea,lbase−1).", "meta": { "altText": "\\mathrm{RRMSE}\\%_{\\mathrm{area},l}=100\\%\\Bigl(\\frac{\\mathrm{RRMSE}_{\\mathrm{area},l}}{\\mathrm{RRMSE}_{\\mathrm{area},l}^{\\mathrm{base}}}-1\\Bigr)." } }, { "type": "Paragraph", "content": [ "When evaluating the forecasts for the entire region as a whole, we use the total RRMSE%, since the areas have very\ndifferent rated power productions. Hence for each\nresolution (", { "type": "MathFragment", "mathLanguage": "mathml", "text": "l", "meta": { "altText": "l" } }, ")" ] }, { "type": "MathBlock", "id": "S5.Ex17", "mathLanguage": "mathml", "text": "RRMSE%ltotal=100%⁢(∑area=115RRMSEarea,l∑area=115RRMSEarea,lbase−1).", "meta": { "altText": "\\mathrm{RRMSE}\\%_{l}^{\\mathrm{total}}=100\\%\\biggl(\\frac{\\sum_{\\mathrm{area}=1}^{15}\\mathrm{RRMSE}_{\\mathrm{area},l}}{\\sum_{\\mathrm{area}=1}^{15}\\mathrm{RRMSE}_{\\mathrm{area},l}^{\\mathrm{base}}}-1\\biggr)." } }, { "type": "Paragraph", "content": [ "Similarly, we will examine the performance for the entire region based on forecast horizon. We previously introduced\nthe ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "RMSEarea,l", "meta": { "altText": "\\mathrm{RMSE}_{\\mathrm{area},l}" } }, " which sums over all prediction horizons. Therefore, we introduce the hourly RMSE as\nwell, which based on the 12 months of results, for a given area, horizon ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "h", "meta": { "altText": "h" } }, " and resolution ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "l", "meta": { "altText": "l" } }, ",\nis given as" ] }, { "type": "MathBlock", "id": "S5.Ex18", "mathLanguage": "mathml", "text": "RRMSEarea,l,h=∑i=1365(yh−24⁢i−y^h−24⁢i)2365.", "meta": { "altText": "\\mathrm{RRMSE}_{\\mathrm{area},l,h}=\\sqrt{\\sum_{i=1}^{365}\\frac{(y_{h-24i}-\\hat{y}_{h-24i})^{2}}{365}}." } }, { "type": "Paragraph", "content": [ "Again, summing across areas,\nwe get" ] }, { "type": "MathBlock", "id": "S5.Ex19", "mathLanguage": "mathml", "text": "RRMSE%l,htotal=100%⁢(∑area=115RRMSEarea,l,h∑area=115RRMSEarea,l,hbase−1).", "meta": { "altText": "\\mathrm{RRMSE}\\%_{l,h}^{\\mathrm{total}}=100\\%\\biggl(\\frac{\\sum_{\\mathrm{area}=1}^{15}\\mathrm{RRMSE}_{\\mathrm{area},l,h}}{\\sum_{\\mathrm{area}=1}^{15}\\mathrm{RRMSE}_{\\mathrm{area},l,h}^{\\mathrm{base}}}-1\\biggr)." } }, { "type": "Paragraph", "id": "S5.p5", "content": [ "For all of these scores, if the RMSE of the reconciled forecast is lower than the RMSE of the base forecast, the\nRRMSE% will be negative, indicating an increase in the accuracy of the forecast. Alternatively, if the accuracy is\ndecreased, the RRMSE% will be positive." ] }, { "type": "Paragraph", "id": "S5.p6", "content": [ "To get a general picture of the performance of the three methods for reconciling wind power forecasts, we examine the\n", { "type": "MathFragment", "mathLanguage": "mathml", "text": "RRMSE%ltotal", "meta": { "altText": "\\mathrm{RRMSE}\\%_{l}^{\\mathrm{total}}" } }, " for all the different temporal resolutions in Table ", { "type": "Cite", "target": "S5-T4", "content": [ "4" ] }, "." ] }, { "type": "Table", "id": "S5-T4", "caption": [ { "type": "Paragraph", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "RRMSE%ltotal", "meta": { "altText": "\\mathrm{RRMSE}\\%_{l}^{\\mathrm{total}}" } }, " of the three methods for each temporal resolution. Greatest reductions are highlighted in bold." ] } ], "rows": [ { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "Temporal resolution [h]" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "MinT-Shrink" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "MinT-Var" ] } ] }, { "type": "TableCell", "content": [ { "type": "Strong", "content": [ "MinT-PVar" ] } ] } ], "rowType": "Header" }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "24" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−19.34", "meta": { "altText": "-19.34" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−19.11", "meta": { "altText": "-19.11" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−19.43", "meta": { "altText": "{-19.43}" } } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "12" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−17.12", "meta": { "altText": "-17.12" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−17.43", "meta": { "altText": "-17.43" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−17.74", "meta": { "altText": "{-17.74}" } } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "8" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−14.96", "meta": { "altText": "-14.96" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−15.13", "meta": { "altText": "-15.13" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−15.38", "meta": { "altText": "{-15.38}" } } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "6" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.71", "meta": { "altText": "-13.71" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.84", "meta": { "altText": "-13.84" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−14.05", "meta": { "altText": "{-14.05}" } } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "4" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.73", "meta": { "altText": "-13.73" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−14.04", "meta": { "altText": "-14.04" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−14.26", "meta": { "altText": "{-14.26}" } } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "3" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.47", "meta": { "altText": "-13.47" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.85", "meta": { "altText": "-13.85" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.99", "meta": { "altText": "{-13.99}" } } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "2" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−12.87", "meta": { "altText": "-12.87" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.29", "meta": { "altText": "-13.29" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−13.46", "meta": { "altText": "{-13.46}" } } ] } ] }, { "type": "TableRow", "cells": [ { "type": "TableCell", "content": [ "1" ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−1.46", "meta": { "altText": "-1.46" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−1.95", "meta": { "altText": "-1.95" } } ] }, { "type": "TableCell", "content": [ { "type": "MathFragment", "mathLanguage": "mathml", "text": "−2.08", "meta": { "altText": "{-2.08}" } } ] } ] } ] }, { "type": "Paragraph", "id": "S5.p7", "content": [ "Table ", { "type": "Cite", "target": "S5-T4", "content": [ "4" ] }, " shows that forecasts on all temporal resolutions are improved, with greater\nimprovements in accuracy for the higher temporal resolutions. As expected, for the 2–24-hour resolutions significant\nimprovements in accuracy are observed when reconciled with the commercial state-of-the-art 1-hour forecasts. However,\nfor the 1-hour forecast resolution, improvements are also seen.\nWith the exception of the 24-hour resolution, the two proposed methods both outperform MinT-Shrink on all resolutions,\nwith MinT-PVar showing the best improvements on all resolutions. The difference between MinT-Var and MinT-PVar is quite\nminor, which tells us that the forecast-dependent variance structure is the main source of the improvements. However,\nadapting the correlation structure also clearly gives an improvement." ] }, { "type": "Paragraph", "id": "S5.p8", "content": [ "Since the base forecasts on the 1-hour resolution are state-of-the-art and used commercially, improvements here are\nsignificantly more valuable for operational purposes than the other resolutions. Therefore, these improvements are\nexamined in further detail in the rest of the section. Figure ", { "type": "Cite", "target": "S5-F6", "content": [ "6" ] }, " shows how the improvements on the\n1-hour resolution are distributed across the 24-hour prediction horizon and across the 15 areas." ] }, { "type": "Figure", "id": "S5-F6", "caption": [ { "type": "Paragraph", "content": [ "Left: ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "RRMSE%l,1⁢htotal", "meta": { "altText": "\\mathrm{RRMSE}\\%_{l,1\\mathrm{h}}^{\\mathrm{total}}" } }, " for each of the three methods across the 24-hour prediction\nperiod. Right: ", { "type": "MathFragment", "mathLanguage": "mathml", "text": "RRMSE%area,1⁢h", "meta": { "altText": "\\mathrm{RRMSE}\\%_{\\mathrm{area},1\\mathrm{h}}" } }, " for each method on each of the 15 areas. The size of the points has\nbeen scaled relative to the rated power production of the area." ] } ], "content": [ { "type": "ImageObject", "contentUrl": "172-f6.png", "mediaType": "image/png", "meta": { "inline": false } } ] }, { "type": "Paragraph", "id": "S5.p9", "content": [ "The three methods that have been tested are able to improve the accuracy of the base forecast across almost all\nhorizons. Only MinT-Shrink falls below the base forecast in accuracy, and only for a single horizon. There seems to be\na pattern to the accuracy improvements for all methods, where the RRMSE drops until the 13-hour horizon and then rises\nuntil the end of the 24-hour cycle." ] }, { "type": "Paragraph", "id": "S5.p10", "content": [ "Looking at how the methods perform for each area, MinT-Var and MinT-PVar once again come out as the best methods. In\nall areas tested, either MinT-Var or MinT-PVar had the best performance. MinT-Var falls behind MinT-Shrink in areas\nno. 4 and 5, while MinT-PVar is consistently better than MinT-Shrink. This highlights the importance of adapting the\ncorrelation to the forecast, as for some areas only adapting the variance is clearly not enough." ] }, { "type": "Paragraph", "id": "S5.p11", "content": [ "Altogether the results show that a reconciliation based on temporal hierarchies can improve even already\nstate-of-the-art commercial wind power forecasts. Although the average improvement of approximately 2% does not seem\nlike much, such improvements are still significant for what are already high-quality forecasts. Further applying the\nvariance estimation techniques proposed in this study results in additional improvements in accuracy, especially when\ncombined with the use of Pearson residuals for correlation estimation." ] }, { "type": "Heading", "id": "S6", "depth": 1, "content": [ "6 Discussion and conclusion" ] }, { "type": "Paragraph", "id": "S6.p1", "content": [ "The investigations performed in this study have shown that, even when using rather simple base forecasts for all the\naggregated levels, commercial state-of-the-art base forecasts can be improved using forecast reconciliation based on\ntemporal hierarchies. This finding is very promising, as integrating the method of forecast reconciliation into\ncommercial forecasts can thus help meet the demand for increasingly accurate forecasts. The additional coherency\nproperty will also be advantageous for the TSOs, as the planning across multiple time horizons becomes simpler and more\ncoherent." ] }, { "type": "Paragraph", "id": "S6.p2", "content": [ "As discussed throughout the study, reconciling forecasts of wind power production poses some challenges in terms of the\nnon-constant variance. The variance will depend greatly on, e.g., the wind speed. Introducing a method which\nadapts the variance structure to the weather forecast thus helps in addressing this challenge. This more rigorous\napproach to reconciliation reflecting the conditional variance dependency on the prediction has resulted in better\naccuracy of the reconciled forecast." ] }, { "type": "Paragraph", "id": "S6.p3", "content": [ "For cases from other fields where data are close to Gaussian, it is more reasonable to assume that the covariance is\nconstant or at least independent of the mean. Hence, we expect less difference between our proposed method and the\nMinT-Shrink in these cases." ] }, { "type": "Paragraph", "id": "S6.p4", "content": [ "There is, however, still much that can be done to build upon the results found here, and as such these results can be\nseen as a step towards an improved methodology for modelling the required covariance structure for reconciling\nforecasts in cases with non-constant variance. Furthermore, this can help build toward a more general understanding of\nthe uncertainty in forecasts when using hierarchies." ] }, { "type": "Paragraph", "id": "S6.p5", "content": [ "Other authors have also developed methods for handling, e.g., bounded data [", { "type": "Cite", "target": "bib-bib18", "content": [ "18" ] }, "]. It would\ntherefore be interesting to examine how this performs compared to our proposed adjustments, or if these methods could\nbe combined in a meaningful way. Due to the limits imposed by the data in this study, further testing on different data\nwould also be beneficial to more confidently show which method is preferable when reconciling forecasts of wind power\nproduction." ] }, { "type": "Paragraph", "id": "S6.p6", "content": [ { "type": "Emphasis", "content": [ "Acknowledgements. " ] }, "\nWe extend a sincere thank you to Krzysztof Burnecki for kindly inviting us to write this article." ] }, { "type": "Paragraph", "id": "S6.p7", "content": [ { "type": "Emphasis", "content": [ "Funding." ] }, " This work is supported by ", { "type": "Emphasis", "content": [ "SEM4Cities" ] }, " (Innovation Fund Denmark, no. 0143-0004), ", { "type": "Emphasis", "content": [ "ARV" ] }, " (H2020-101036723),\n", { "type": "Emphasis", "content": [ "ELEXIA" ] }, " (EU HE-101075656), ", { "type": "Emphasis", "content": [ "IEA Wind Task 51" ] }, " (EUDP nr. 134-22015), and finally the projects ", { "type": "Emphasis", "content": [ "PtX,\nSector Coupling and LCA" ] }, " and ", { "type": "Emphasis", "content": [ "DynFlex" ] }, ", which both are part of the Danish Mission Green Fuel portfolio\nof projects." ] }, { "type": "Paragraph", "id": "authorinfo", "content": [ "\nMikkel Lindstrøm Sørensen has an MSc in mathematical modelling and computing from the Technical University of Denmark.\nHe is currently a PhD student working on renewable energy forecasting and hierarchical forecast reconciliation, under\nthe supervision of Jan Møller and Henrik Madsen.\n", { "type": "Link", "target": "mailto:mliso@dtu.dk", "content": [ "mliso@dtu.dk" ] }, "\nJan Kloppenborg Møller is an associate professor in stochastic dynamical systems at the Technical\nUniversity of Denmark. He holds an MSc of applied mathematics (2006) and a PhD in engineering from\nDTU (2011). His research is concentrated on modelling and forecasting of (continuous- or discrete-time) stochastic\ndynamical systems.\n", { "type": "Link", "target": "mailto:jkmo@dtu.dk", "content": [ "jkmo@dtu.dk" ] }, "\nHenrik Madsen received his PhD from the Technical University of Denmark (DTU) in 1986. Since then he has been employed\nby DTU, first as assistant professor, then associate professor from 1989, and he has been a full professor in\nstochastic systems since 1999. In June 2016, he was appointed as Knight of the Order of Dannebrog by Her Majesty the Queen\nof Denmark, and he was appointed as Doctor h.c. at Lund University in June 2017. His research areas are primarily\nmodelling, forecasting and control, with a special focus on applications related to dynamics and\nflexibility in energy systems.\n", { "type": "Link", "target": "mailto:hmad@dtu.dk", "content": [ "hmad@dtu.dk" ] } ] } ] }