Bayesian framework for power network planning under uncertainty

Effective transmission expansion planning is necessary to ensure a power system can satisfy all demand both reliably and economically. However, at the time reinforcement decisions are made many elements of the future system background are uncertain, such as demand level, type and location of installed generators, and plant availability statistics. Making decisions which account for such uncertainties is presentlyusuallydonebyconsideringasmallsetofplausiblescenarios,andtheresultinglimitedcoverage of parameter space limits confidence that the resulting decision will be a good one with respect to the real world. This paper presents a methodology which uses statistical emulators to quantify uncertainty in mathematical model outputs for all points at which it has not been evaluated, and hence to control properly uncertainties in the decision process arising from the finite size set of scenarios. The result is a generally applicable approach to network planning under uncertainty, including decision makers’ risk preferences, which scales well with problem size. The approach is demonstrated on a Great Britain test problem, which replicates key features of the model the Transmission Owners use for practical strategic planning studies.


Introduction
Transmission system planning is key to ensuring that a system's future demand can be met both economically and with an appropriately low adequacy risk. Historically, systems were planned to maintain continuity of supply with any α components on outage, α commonly being 1; this is known as an N-α planning criterion [1][2][3][4]. More recently, a number of authors (e.g. [1,5,6]) have suggested that the likelihood of events occurring must also be considered, because it would typically be deemed undesirable to make a large investment to protect against a very unlikely event, or one whose consequences are very small. Further, in an economic setting, as will be presented in this paper, it would be similarly undesirable to make a decision that alleviates costs in a handful of negligible scenarios, whilst resulting in a significant over-investment in any scenario that could plausibly occur. Therefore, it would be more relevant to consider problems which not only satisfy the N-α criterion, but also make decisions to optimise some secondary economic/welfare function.
Such secondary criteria may include: minimising the sum of construction costs, operation costs and stand-by costs [7]; minimising the present worth of investment costs and expected outage and production costs [8]; or minimising the sum of investment costs and operation costs [9]. [10] demonstrates a holistic view which bases decisions on a weighted average of multiple criteria across multiple stakeholders. A cost-benefit analysis (CBA) approach can also be taken [11] (to optimise incremental costs plus energy replacement costs) and [12] (to optimise the ratio of benefits a reinforcement brings against the cost of the reinforcement). [2] considers basing decisions on investment sensitivities, i.e. ratio of improving the objective function (e.g. satisfying of the supply-demand curves or power exchange deviations) against investment cost. An approach using co-operative game theory is used by [13], involving creating coalitions between multiple agents to both increase the benefits for each agent and reduce the number of lines needed to expand the network. [14] considers the marginal profits of reinforcements made against the marginal cost to reinforce whilst accounting for the intermittent generation associated with renewable resources. [15] considers reinforcement planning to minimise adjusted delivered cost (delivered cost of a resource to a load zone which considers bus-bar and transmission (investment, operation and line losses) costs and also adjusts for key market value factors (integration costs, avoided resource adequacy costs and avoided http time-of-delivery costs)). Further, it considers how cost estimates and resulting decisions change as input factors are varied (such factors are the proportion of renewable generation on the system and the associated costs with renewable generation, such as tax incentives).
When making planning decisions it is important to account for the uncertainty that exists in background against which the system is planned. For example, if demand is higher than expected there may be insufficient transmission capability to satisfy it, or if it is lower some assets may with hindsight be deemed stranded. This is commonly handled by making the decision considering a limited number of scenarios, e.g. [16] which formulates a mixed integer programming model to minimise the total expected cost of generator operation and investment and transmission investment across several scenarios. The expected value of perfect information, and the expected cost of ignoring uncertainty, are calculated to quantify the potential benefits of better information and costs of not considering uncertainty in system background. [17] maximises generation company and transmission company expected profits, subject to a reliability check from the independent system operator which accounts for the uncertainty that exists in future load growth and outages. [18] seeks to minimise expansion investment whilst maximising system reliability and security for each of several scenarios which vary system load and additional capacity installations. Initially, each scenario is optimised separately, then the expansion with the lowest adaptation cost (i.e. the cost for additional reinforcement if the initial scenario is incorrect) is selected. [19] seeks to make a decision which accounts for the uncertainty that exists due to the varying availability of hydro power. This decision minimises the cost of investment, plus a weighted average of expected costs due to load shedding across multiple scenarios.
However, a common feature to all these methodologies is that uncertainty is handled by considering a finite (and usually small) set of possible scenarios. Particularly where the uncertainty in planning is high dimensional, such a finite set of scenarios may give sparse coverage of the space of model inputs, and hence an optimal decision within the model which is strongly dependent on the particular choice of scenarios rather than the underlying uncertainty in model inputs.
The key contribution of this paper is to demonstrate how the use of a statistical emulator, which quantifies the uncertainty in the decision support model outputs for inputs where it has not been evaluated, and hence makes possible a systematic broad search of the decision space with relatively few model evaluations. The exemplar used to demonstrate this emulation approach is closely based on the probability model used by National Grid in practical Great Britain planning studies, and is described in Section 2. It performs a cost-benefit analysis (CBA) between network capital costs, and expected constraint costs arising from redispatch of generation due to finite network capacity.
Using the first emulator covering the whole decision space, the range of the decision space which cannot yet be rejected as sub-optimal is still quite large (due to large uncertainty in the model output at points far from model evaluations). However by performing further model runs in the area of decision space which remains of interest, it is possible to fit an improved emulator there, and hence narrow down the search space further. This process of narrowing down the area of search space considered, through multiple waves of emulation, makes the approach scalable to much larger problems than the exemplar considered here.
The emulation method is described in detail in Section 3construction of the emulator is demonstrated on the Great Britain planning example in Section 4, and its application in network reinforcement decisions is demonstrated in Section 5. The sensitivity of decisions to precise quantification of prior knowledge, and to attitude to risk, are explored in Sections 6 and 7 respectively.

Problem specification
This paper considers the problem of transmission expansion planning from an economic perspective, in which decisions are taken based on a cost-benefit analysis between the capital costs of network expansion and future constraint costs (the costs of finite network capacity constraining the generation schedule). The additional network capacity is a decision variable, and the capital costs of a given upgrade are usually known to a high degree of accuracy. In contrast, the assessment of future constraint costs must be carried out under considerable uncertainty regarding system background (plan availability statistics, demand level, generation prices etc.).

Constraint cost model and Great Britain data
The forward energy market in GB is not locational, i.e. it is run as if there were infinite network capacity. After 'gate closure' one hour ahead of real time the System Operator (National Grid) must then adjust the generation schedule in order to balance aggregate supply and demand, and to ensure a feasible network flow. In GB the big picture of network flow patterns is from north to south, so for strategic planning it is appropriate to model the transmission network using a tree-like network model consisting of regions separated by boundaries of finite capacity, as illustrated in Fig. 1. This paper will use as an example reinforcement of the B6 and B7a boundaries, which lie around Northern England; this has been one of the most important reinforcement issues in recent years [20].
Given fixed input data 1 the out-turn constraint costs in a future year may be calculated as where t indexes the periods into which the year is subdivided; z t is a vector of random variables representing available capacities of each generating technology in each region in period t and the vector s t is the input demands in each region in period t. The function g evaluates the constraint costs in a period given demands and generation availabilities, by formulating the redispatch as a linear program. Following National Grid's practice, we evaluate the expected constraint costs E[C ] in the future year; due to the standard result for the expected value of a sum of random variables, to evaluate E[C ] it is not necessary to consider serial associations between generation availabilities in different periods. This model output is evaluated by Monte Carlo simulation, using importance sampling [21] to improve convergence by sampling more densely demand levels which contribute most strongly to the result. Data are taken from the Electricity Scenarios Illustrator (ELSI) published by National Grid on their consultation and engagement website [22].

Structure of uncertainty problem
The constraint cost model may be thought of as a function y c = f c (x) (2) where y c represents the expected constraint costs, and x the inputs to the calculation. As f is evaluated by simulation, in principle y is not known with certainty even for fixed x; however in practice for the example used here it is possible to use a sufficiently large Monte Carlo sample that y may be assumed to be known exactly.
In general, the model inputs x are not all known precisely. f c (x) may thus be thought of as f c (v, a, d), where v represents the inputs in which uncertainty is modelled explicitly; a represents the inputs which are either known with certainty or in which uncertainties are not of interest, so their values are treated as fixed as if they are known precisely; and d represents the decision variables.
When evaluating expected constraint costs, it is important to consider uncertainty in input parameters such as demand level and plant availability statistics. This paper will consider how cost estimation and reinforcement decisions are affected by the uncertainties expressed in prior beliefs about nuclear and CCGT availability probabilities as well as peak demand level, parameters to which expected constraint costs have been found to be particularly sensitive. These would thus make up v. a would comprise of installed generation capacities, transfer capabilities on boundaries not considered in the reinforcement decision problem, and other inputs such as energy prices and other generation availabilities. d would be the magnitude of B6 and B7a reinforcement.
In this paper we wish to minimise the total costs of a reinforcement decision, so the costs of concern are i.e. the constraint costs that arise given a particular reinforcement decision plus the cost of the reinforcement itself (f r (d)).

The emulation process
In the problem presented, there is uncertainty in the input values of the inputs x. However, the function, f , is far too expensive to evaluate at every set of input parameters desired. Therefore, it is necessary to evaluate f a small number of times for particular values of x, and approximate f everywhere else as carefully as possible by some alternative function,f .f should be much less computationally demanding to evaluate in order to allow for efficient estimation for input values not simulated. The full emulation model fitted also details the uncertainty that exists in the response when approximating f byf .
Suppose there are k variables with uncertainty of interest, labelled v 1 , . . . , v k and m decision variables, d 1 , . . . , d m The goal is to fit a model to the value of y based on these k + m variables.
As this process assumes all other variables are considered known, The emulator used in this paper will use polynomial regression as the basis for modelling y. This can be thought of as follows i.e. the dot product of a coefficient vector with the polynomial form of the predictor variables with error term ε(v 1 , . . . , v k , d 1 , . . . , d m ).
The estimates of the model parameters,β, have an associated covariance matrix, c(β,β). The model parameters are thus multivariate normally distributed with meanβ and covariance matrix c(β,β). Whilst this methodology applies very generally, the terms included in the regression model should be carefully selected based on the specifics of the example considered. This is considered in Section 3.5.
The fit of the model can be improved further by using a Gaussian process model to smooth the residuals, This is important as the polynomial model will aim to give a good fit over the entire range of the variables, whereas the Gaussian process can model local behaviour much more accurately, as described below. The idea of a Gaussian process is to smooth the residuals in order for the model to agree with the training data and provide an accurate estimate of response (and variance of the response) for values where training data is unavailable. In this sense, the Gaussian process model can be thought of as a smooth interpolator of the residuals [23]. In simple terms, the idea of the Gaussian process is if a prediction at point x p is to be taken, training data closer to x p carries more relevant information and should be given more weight when making a prediction.
Gaussian Process models fitted to a set of training data, in give a mean-variance relationship (i.e. expected response and uncertainty in that estimate) for the response as a function of each possible input value [23,24]. [23] gives an accessible introduction to Gaussian processes models. The exemplar model is a simple function of a single variable (specifically x + 3 sin( x 2 )), acting in place of a complicated simulator for demonstration purposes. The function is evaluated exactly a small number of times (training runs) and it is demonstrated how a smooth interpolation (the Gaussian process) of these training runs acts as a good approximation to the true function, as well as illustrating the approximation error at inputs where the model has not been evaluated. This tutorial reference also demonstrates how the Gaussian process emulator can improved by using additional training runs, and how the parameters of the Gaussian process itself affect the quality of the resulting approximation. [25] is a more technical reference on fitting Gaussian process models and their use in making predictions. In order to understand how to fit and make predictions from a Gaussian process, the correlation function, κ, must first be defined. For two vectors of predictor variables, , the form of the correlation function for given correlation parameters γ = (γ 1 , . . . , γ k+m ) considered in this paper is In practice the correlation between 2 matrices is required. Suppose X 1 and X 2 are two matrices where each row of each matrix is a set of predictor variables. Further, suppose X 1 has I rows, x 1,i is the ith row of X 1 , X 2 has J rows and x 2,j is the jth row of X 2 . The correlation matrix between these two matrices would be K X 1 ,X 2 , where the jth entry of the ith row of To fit a Gaussian process, the design matrix, X d , of the training runs is required, as well as a vector of calculated responses, ε. The ith row of X d describes the input values used for the ith training run, and the ith entry of ε is the corresponding residual of that training run after the polynomial regression model has been fitted. Suppose prediction is desired at a new set of data-points. A matrix X * is created, where each row details the input values for which a prediction is required. Note, X * can be a single set of input values and will be treated as if it were a matrix with a single row. The vector for estimated response of the Gaussian process for the input values of X * is [25,26] with the covariance matrix of these predictions The parameters γ 1 , . . . , γ n in Eq. (5) control how much weight is given to each training run when making predictions. The larger the values of these parameters, the faster the relative weight given to a training runs decreases with distance from the point where a prediction is required. In simple terms, larger values of these parameters will give a function more sensitive to the local behaviour of training runs than smaller values. A method utilising cross-validation should be used when selecting the Gaussian process model parameters [25], such as the leave one out estimate which will be described in Section 4.2. This is because such methods estimate how well a given model will perform when predicting a response for a new set of input values. This is particularly important for Gaussian Process models, as any choice of model parameters will interpolate the training data used to fit them, so traditional model selection criteria, such as minimising residual square error, would be meaningless.
For values where training data was available, the mean of the response will be the residual of the polynomial model with 0 variance. This means that the training data is interpolated perfectly with zero variance at the training data input. cov(ε * ) is greater than zero away from values where training data input was available, and generally larger the further a set of input is from a point which had training data available.
In this paper the polynomial portion of the emulator will be noted byf 1 and the Gaussian process model fitted to its residuals will be noted byf 2 . Thus, the emulator model for given inputs is the sum of these two, denotedg, such that That is to sayg is our emulator which approximates the simulator, and the emulator estimate is acquired by evaluating the sum of the fitted polynomial regression model and corresponding Gaussian process model for given input values.

Cost estimation from emulation
The expectation of a function of a continuous random variable is calculated as the integral of the product of the function and the probability density function (PDF) of the random variable. Applying this to emulation, if the PDF of the uncertain variables (i.e. v 1 , . . . , v k ) was known to be p(v 1 , . . . , v k ), the estimate of expected costs for a given decision, d * In real problems, such as the constraint cost problem presented, the PDF of the uncertain variables is often not known exactly. In particular, in cases such as this, assessment of uncertainty in the planning background relies to a large extent on expert judgement [27,28], and in a Bayesian formulation the prior judgement of uncertainty should be expressed as a PDF of the variables v 1 , . . . , v k . The expected costs under uncertainty can then be estimated as in Eq. (9). Examples of prior beliefs for the constraint cost problem are given in Section 4.1.
When making a decision it is desirable to identify the value of the decision variable which minimises/maximises (as appropriate) Eq. (9). This means that the integral in Eq. (9) must be evaluated multiple times. When d 1 , . . . , d m are discrete and a small number of possible decisions exist this is straightforward, as one can simply evaluate explicitly all possible decisions and select the best one.
When d 1 , . . . , d m are continuous or can take a large number of discrete values, numerical methods may be required to identify the optimal decision, as we will demonstrate in the example.
As mentioned in Section 3.1, there also exists uncertainty in the fitted emulation model itself. Recall the polynomial portion of the model,f 1 , has model parameters β. There is uncertainty in the estimates of the values of these parameters and our model of their true value is said to be a multivariate normal distribution with meañ β with covariance matrix c(β,β). A new set of model parameters can then be randomly drawn from this multivariate normal distribution, call the corresponding polynomial model using such a random drawing of model parametersf 1,r . Corresponding residuals, ε r , can then be calculated for usingf 1,r instead off 1 when predicting. A new Gaussian process model,f 2,r , can then be fitted as in Section 3.1, using the error residuals of this randomly drawn polynomial model, ε r , in place of the error residuals of the original polynomial model. A random variation of the emulator model could thus be the sum of the polynomial model with a random draw of parameters and the corresponding Gaussian process fitted to the resulting residuals such that: However, such variations in the model parameters can have a large effect on the fitted model. Further, usingg r instead ofg in Eq. (9) could also have a large effect on the estimated costs. By carefully considering how the cost estimates vary as the model parameters are varied, credible bands can be formed which give a range in which we would expect the actual costs to lie for a given level of confidence.
To do this, consider creating a randomly drawn model, g r , as described above. For a given decision, d * 1 , . . . , d * m , an expected response using the randomly drawn model can be calculated by integrating over uncertainties, as in Eq. (9) such that:

Latin hypercubes
Training runs to be used in the emulation process were acquired using Latin hypercube (LHC) sampling [29]. When taking a Latin hypercube sample, first the size of the sample, n, must be specified. Then each axis of the hypercube formed by the variables of interest (i.e. v 1 , . . . , v k , d 1 , . . . , d m ) is divided into n intervals. A sample of size n is then taken which varies the values of the k + m variables such that for each variable exactly one of the samples has a corresponding value in each of the n intervals of that variable's axis. [29] includes illustrations to demonstrate this process.
The Latin hypercube sample gives the input values to be used when simulating training runs, which in turn are used to fit the emulator. An appropriate number of training runs must be selected and is considered in Section 4.2 via cross validation.
Latin hypercubes are very advantageous in comparison to grids (i.e. evenly distributed points over a given range). Firstly, a more dense coverage of the region sampled is given, which should allow models to give a better fit. Secondly, significantly fewer points can be used when making decisions.
In the example presented in Section 4 it will be shown only a 300 point LHC sample is necessary, whereas a 5 by 5 by 5 by 5 by 5 by 5 grid would use over 15,000 points and the resulting emulator would not give a particularly good fit. This saves a considerable amount of work in the simulation of training runs. This benefit would be even greater in higher dimensions, where grid size increases exponentially with dimension, whereas LHC samples do not.

Reinforcement decision making
The objective is to make a reinforcement decision which minimises the total costs in the system (the sum of constraint costs plus the costs of reinforcement). This is achieved by taking waves of observations, where in each wave the aim is to eliminate potential decisions if there is sufficient evidence to suggests they are not optimal in the following manner: 1. Set an initial range of values to be considered for the decision variables. 2. Simulate training runs based off a Latin hypercube sample using the current limits on the decision variables. 3. Fit an emulation model to the latest simulated sample and use it to eliminate decisions (details below). 4. Repeat until uncertainty about optimal decision reaches an acceptable level. The optimal decision is the decision which minimises Eq. (9).
To do this, a LHC sample will be taken over a specified initial range of reinforcements to consider (in this example this is between 0 and 4000 MW for both B6 and B7a). This LHC sample gives the values of input variables to be used for the training data. Simulations are then taken using these input values to give the training data. Then the emulation process of Section 3.1 uses this training data to create an emulator model of how input affects output of the simulator. against that decision being optimal, as the lower bound of the estimate is not better than the best upper bound. That decision can thus be eliminated from the decision space.
A new set of training runs is taken where the values of the decision variables are only those not eliminated. A new emulator model is fitted to this new set of training runs, and credible intervals are estimated for decisions which were not previously ruled out. This allows for further bad decisions to be identified and eliminated. The process continues iteratively until the uncertainty about the optimal decision reaches an acceptable level.
The use of a wave process is actually very efficient. Initially, global behaviour is modelled, so an area where the optimal decision is likely to lie can be identified. Subsequent waves can then model local behaviour much more accurately meaning an accurate decision can be made. Further, a much better model (and therefore better decision) can be obtained from several waves of relatively small number of observations than a single wave of a relatively large number of observations. Reinforcement decisions are typically considered over multiple years. For simplicity, in this paper details have not been given for the emulation process over multiple years. However, the simulation process presented naturally extends to multiple years, by detailing the same values of x for each year to be simulated and calculating y as the sum of all expected constraint costs over all periods of all years plus the cost of any reinforcements built.

The particular example considered
The emulation example presented in this paper considers uncertainty in 3 input values: nuclear availability probability, CCGT availability probability and peak demand value. These are considered to be variables v 1 , v 2 and v 3 respectively. This means other input values will be kept constant. Two decision variables are considered: B6 and B7a reinforcement magnitudes. These are considered to be variables d 1 and d 2 respectively. These variables were chosen as a result of preliminary experimentation into which variables the cost estimates appear to depend most sensitively upon. Further the B6 and B7a boundaries are adjacent boundaries to the North of England/South of Scotland, which are important in the time frame of the example as a lot of wind generation is scheduled to be built in Scotland.
An emulator will be fitted to model how the total costs in the system (constraint costs plus the costs of any reinforcement made) are affected by the values taken by these 5 parameters. The polynomial portion of the model includes first and second order terms of all 5 variables. Interaction terms (pairwise and threeway) are included between v 1 , d 1 and d 2 as well as interaction terms (pairwise and three-way) of v 2 1 , d 2 1 and d 2 2 . The residuals of this polynomial are used to fit a Gaussian process depending on variables v 1 , v 2 , v 3 , d 1 and d 2 as described in Section 3.1. The sum of these two models forms our emulator, as in Eq. (8). Table 1 summarises the parameters included in the linear model, along with corresponding value of model coefficient (the particular β from Eq. (4)). The model coefficients of the Gaussian process model, γ 1 , . . . , γ 5 in Eq. (5), are 20, 270, 1.5 × 1 −7 , 1.9 × 10 −5 and 2.6×10 −6 respectively. These are all given to 2 significant figures.
Cross-validation, detail given for sample size selection in Section 4.2, was used to select the terms to be included in 1.9 × 10 2 the polynomial regression model. As B6 and B7a are adjacent boundaries it was natural to expect an interaction between the two. Further, a high merit capacity such as nuclear will incur large constraint costs if it is unable to be used, with these costs being reduced heavily as transmission capacity is increased, explaining its inclusion in the interactions. The R 2 statistic of the fitted model can also be useful for model selection, in particular for being a good quick reference for how high degree a polynomial model seems necessary. Analysis of variance (anova) is also a useful tool to ensure model parameters do indeed add explanatory power to the model. The use of prior beliefs and the methodology of Section 3.2 can then be used to estimate costs under uncertainty as the decision is varied, with the goal being to identify the decision which minimises the expected costs in the system. While this example has fairly low dimensional input space, the methodology used for managing uncertainty scales well to a much higher dimensional input and decision space, for instance [30] where up to 17 uncertain parameters are considered.

Uncertain variables and prior beliefs
The prior beliefs used for illustration in this paper are based on the expert judgement of Paul Plumptre (PP), formerly of National Grid [31]. For the three uncertain variables, a range is given for values the variables could feasibly take. For Nuclear availability probability, the given range is 0.55 and 0.85 and for CCGT the given range is 0.8-0.95. PP also states an appropriate level of uncertainty in peak demand level is 1% for each future year used. Therefore, the peak demand level will be considered to be multiplied by a value between 0.95 and 1.05. A uniform distribution is used to represent prior beliefs across each of these ranges.
Year 6 data from [22] is used as the basis of the other aspects of the power system for the exemplar. This means input parameters for the power system which are treated as fixed (such as the load duration curves, other availability probabilities, installed generation capacities etc.) will be taken as the values given by [22]. This is intended to represent a 2016 power system scenario. The exception is the boundary capabilities of B6 and B7a will be set to year 1 levels. This is as if a reinforcement decision for 2016 power system is being made based on the estimates of the 2016 system from information available in 2010. Other boundaries will be treated as if they have infinite capacity so only costs arising at B6 and B7a are considered. Reinforcements between 0 and 4000 MW on each boundary are considered.

Training run selection
Cross-validation is used to assess how many training runs are necessary to fit an emulator. This is a criterion based on the ''leave where n is the number of points in the sample. Table 2 displays how the cross-validation values vary with the number of points used in the sample. It is desirable to minimise the value of the cross-validation, as this indicates a better prediction for input values not simulated.
It can be seen that the value of cross-validation gradually decreases with sample size. There is a noticeable improvement in increasing from a 100 to a 200 point Latin hypercube sample. There is a smaller, but noticeable improvement increasing from a 200 to a 300 point Latin hypercube sample and further improvement from 300 point to 400 point sample. However, there does not appear to any gain in using a sample size larger than this.
The R 2 value of the polynomial proportion of the model is also given. All values are above 0.97 which indicates an excellent fit is given, even without the inclusion of a Gaussian process model.

Graphical illustration of uncertainty
Fig. 2 displays how the shape of the surface of constraint costs varies with level of peak demand when varying nuclear and CCGT availability probabilities. These are exact values from simulation, not emulation models, and assume no reinforcement has been made. There is not a great deal of difference between the plots, indicating that the demand level has very little effect on results. It is also seen that the nuclear availability probability appears to have the largest effect on estimates, which suggest cost estimates will be most sensitive to prior beliefs about it.
Cross-validation was again used to assess which model terms are necessary to include. Although the graphs were relatively flat, it appeared necessary to include quadratic terms. Further, interaction terms were necessary between the two decision variables and nuclear availability probability, though interactions involving other variables did not seem necessary.   both boundaries reinforcement are above 2500 MW. This indicates all decisions near the optimal solution give very similar results. Section 3.2 indicates that all estimates come with upper and lower bounds. Further, Section 3.4 indicates how these bounds can be used to eliminate bad decisions, allowing a more accurate model to be fitted over a smaller range of decisions. This allows for a better, more confident decision to made. The improvement in the model is explored further in Section 5.2. Fig. 4(a) shows which decisions seemed to merit further consideration in the second wave (i.e. based on the methodology of Section 3.4 there was no conclusive evidence of them being worse than the apparent optimal). Values of the B6 reinforcement magnitude range from 2300 to 4000 MW and values of the B7a reinforcement range from 1900 to 4000 MW. This is still a quite large region, though it has eliminated 82.7% of the decisions considered in the first wave.

First wave and further waves
Note, when sampling training runs to use in the next wave the shape of Fig. 4 is preserved. First, the decision which gives the minimum upper bound in the current wave is identified as , as in Section 3.4. Then, a set of values for the decision variables are randomly drawn from a uniform distribution. A lower bound for the estimate using these decision is evaluated using methodology from Section 3.2. If this lower bound is less than ) the values of the decision variables are used for training runs for the next wave. If not, the decision is rejected and new values are selected for the decision variables. This is repeated iteratively until the specified number of training runs have been acquired. Values for the variables whose values are uncertain are sampled via a Latin hypercube as described in Section 3.3. A new emulator model is then fitted to the new set of training runs. This allows more decisions to be eliminated by considering the bounds of the estimates, as in Section 3.4. Fig. 4(b) displays the reinforcement combinations that would merit further consideration in a third wave of reinforcements. Values of reinforcement magnitude of the B6 boundary now range from 2800 to 3400 MW and values of B7a reinforcement magnitude ranges from 2300 to 3100 MW. If sufficient resources are available, the wave process can be continued, eliminating even more decisions.
Alternatively, after sufficient decisions have been eliminated a final decision could be taken to minimise costs for the current wave. The reason to make such a decision is the uncertainty in a model will not always allow for all but a single decision to be eliminated but such precise decisions are useful for real life implementation of reinforcement. Note, this would not be a definitive optimum, as uncertainty is always present in the model and consequently any decision not eliminated in this wave could also potentially be optimal.
If such a decision was taken in the first wave this would result in 3370 MW B6 reinforcement and 2830 MW B7a reinforcement. However, such a decision in the third wave would result in 3060 MW B6 reinforcement and a 2890 MW B7a reinforcement, somewhat different to the decision that would be made in the first wave indicating the improvement in the model. The rest of this paper will consider taking a final decision to minimise costs in the third wave.  Fig. 5(a). It can be seen that the two models are very different, with costs varying much more smoothly in the third wave. This is an improvement due to the fact that the third wave model is fitted over a much smaller range of decisions, allowing it to be much more accurate over that range. Fig. 6(a) illustrates the construction of credible intervals for the first wave as the B6 reinforcement is varied. In this graph B7a reinforcement is fixed to the apparent optimum of 2890 MW. 200 grey curves are shown, each of which gives the resulting estimate of costs from random variations of the fitted model, using the methodology of Section 3.2. Credible intervals for the model are shown with dashed black lines. It can be seen that these estimates capture the majority of the variation of the estimates. Further, the expected model, shown in solid black, is contained within these limits. Fig. 6(b) illustrates the improvement in the credible interval of the estimate between the third and first wave. The first thing to note is the model from the third wave is much smoother than that of the first, indicating the improved fit from fitting over a smaller range of decisions. Additionally, the credible intervals are much narrower for the third wave than the first. This indicates the increased certainty within the estimate. This increased certainty allows for more decisions to be eliminated, or if a decision is to be taken it allows for a narrower range to be given that could possibly contain the optimal.

Final wave results
The power of the presented methodology of emulation through several waves is demonstrated by the fact that a much better model is achieved by taking 3 separate waves of 300 observations, with the final one being over a much narrower range, than attempting to fit one model to a single set of 900 training runs. In this final wave the fitted emulator model is a very accurate approximation of the simulator, whilst being much less expensive to evaluate.

Sensitivity to reinforcement cost
So far, it has been assumed that it costs £1000 per km per MW to reinforce, with 100 km boundary thickness assumed for both boundaries. However, estimates of this cost vary, even within the  same organisation, with [32] stating the value to be £750 per km per MW, and [31] gives the assumed £1000 per km per MW. Table 3 displays how the decisions giving the lowest total costs varies with the input cost for reinforcement. If the cost to reinforce is increased 50% to £1500 per km per MW, the total reinforcement is reduced by 1060 MW. Conversely, if the cost was reduced by 50% to £500 per km per MW, the total reinforcement is increased by 1160 MW.
There also appears to be a smooth general trend of 100 MW reinforcement less being made on each boundary for each extra £100 per MW per km it costs to reinforce. This is to be expected, as there is more economical benefit to building a reinforcement when the reinforcement itself costs less to build.

Varying the prior beliefs
A single set of prior beliefs have been assumed so far in this article. These beliefs were based on the expert judgement of Paul Plumptre, formerly of National Grid. However, the limits given could be expanded or contracted. Varying these two factors has an effect on both the cost estimates and the resulting decisions. This section presents results when varying the intervals given for the prior beliefs. The prior beliefs given in Section 4.1 will be referred to as the PP prior. Fig. 7 displays how the estimated yearly constraint costs vary when the range of the variables is varied from the values given in the PP prior. This graph assumed no reinforcement was made. One constraint put in place is the availability probabilities were capped at 1, as the probability of an event cannot exceed 1. For example, the original prior gave a window for CCGT availability probability as 0.8-0.95. Expanding the window by a factor of 2 would give the width 0.725-1, not 0.725-1.025.
As can be seen, the prior beliefs can have an effect on estimated costs. In particular, costs appear to rise as the width of the uncertainties increases with an increase from £2.35×10 8 to £2.6× 10 8 when uncertainty is increased in all variables. This indicates that if the prior belief intervals are too narrow, there is a risk of underestimating costs and making a suboptimal decision. Table 4 shows how the optimal decision varies with the prior beliefs used. Using a narrower set of beliefs has a small effect, with the B6 boundary being decreased by 40 MW and the B7a boundary being decreased by 210 MW. However, using the wider set of prior beliefs has no effect on the decision made. This indicates that in  this particular example there is little risk of making a bad decision if the prior belief intervals are too narrow/wide. However, this is a conclusion that we are unable to reach once such an analysis has been performed. Further, we earlier observed sensitivity when considering cost to reinforce in Section 6.1 and sensitivity to attitude to risk will be demonstrated in Section 7, indicating the importance to carry out such analysis.

Attitude to risk
Results and methodologies presented so far have assumed a risk neutral attitude. However, it is often desirable to take an attitude to risk into account when making planning decisions, as investors commonly weight negative outcomes more heavily than positive benefits, and in particular may have a particular aversion to very extreme adverse scenarios.

Attitude to risk: methodology
The methodology in this paper can be adapted to allow decisions to be made which take into account an attitude to risk. In order to do this, first the optimal values of the decision variables, A modelg b is then fitted such that: g b is the expected benefit relative to the risk neutral optimal when uncertain variables take values v 1 , . . . , v k and decision d 1 , . . . , d m is made. When considering an attitude to risk it is wise to work with g b and not directly with the costs. This is because in the presented example extremely large costs (hundreds of millions) occur no matter which decision is made. The variation in costs as the values of the uncertain variables are varied is quite small in comparison (tens of millions). This means a utility function would not make a significant change to the decision unless an extreme attitude to risk was taken, despite the variation in costs being very significant in real terms. Working with g b allows for the relative benefits and losses as the uncertain variables are varied when making a particular decision to be adequately assessed. As a result, the utility function can much more accurately reflect the preferences of a real life decision maker.
For a utility function, u, which describes an attitude to risk, the expected utility of a decision can then be calculated as: The objective is to maximise Eq. (14). It is necessary to carry out a separate wave process for each attitude to risk. This is because risk averse attitudes would consider decisions of larger magnitude than a risk neutral attitude, so it is better to work with several models fitted more accurately over a smaller range than one model fitted less accurately over a large range.
In practice, real world decision makers are rarely risk prone. Fig. 8 illustrates this utility function.
Under risk averse conditions, this utility function will harshly penalise a scenario for doing worse than the risk neutral optimal, thus placing in the decision process great emphasis on avoiding scenarios where extremely large losses occur. The utility function also rewards improvement on the risk neutral optimal, however this reward is relatively smaller than the penalty for doing worse than the risk-neutral optimum. Tables 5 and 6 shows how the decision made varies as the attitude to risk is varied when assuming a cost of £1000 and  This example has illustrated that attitude to risk can be an important factor to consider when making a decision. Even in the less sensitive £1000 per MW per km case it was still necessary to do this analysis, as this is the only way to know that the result would not be very sensitive to attitude to risk.

Conclusion
This paper has presented a statistical methodology for taking power system upgrade decisions under uncertainty over the future planning background. Investment decisions are based on a cost benefit analysis between the capital cost of the upgrade and the constraint costs arising from finite network capacity restricting the generation schedule, and a computer model is used to estimate expected constraint costs for a given background and investment decision. Uncertainty in inputs is represented as a probability distribution in a Bayesian framework.
Because of the substantial computational expense of evaluations, it is only possible to evaluate the decision support model a limited number of times, and as a consequence the value of the model output is uncertain almost everywhere in the decision/input space. Key to the approach presented is therefore to use a statistical emulator to quantify this uncertainty at points where the model has not been evaluated. This provides a systematic means of considering the entire decision space even with relatively few model evaluations, and the decision space may efficiently be narrowed down by successively rejecting decisions which are very unlikely to be optimal based on evidence from the current emulator, and constructing an improved emulator through additional model evaluations in the region of decision space which remains of interest.
The approach is demonstrated on a Great Britain test problem which replicates National Grid's constraint cost estimation approach. The effect on investment decisions of reinforcement cost, degree of uncertainty about system background, and attitude to risk is explored -larger reinforcements are optimal given greater uncertainty in model inputs, and given a greater aversion to risk. It is also notable that at higher reinforcement costs results are more sensitive to attitude to risk.
The key advance over previous approaches to related problems is the way in which the emulation approach allows systematic exploration of the whole decision space and space of uncertain inputs, without concern that the optimal solution within the model might depend on the precise choice of discrete scenarios used to represent uncertainty. For some combinations of inputs and attitude to risk, this uncertainty analysis gives an investment plan which is very different from that which would be obtained with a single scenario. Where the difference arising from the uncertainty analysis is smaller (e.g. at lower reinforcement cost), it is only possible to know this by doing the uncertainty analysis to link the model to the real system, and the detailed uncertainty analysis gives greater confidence that even small savings seen in the model will be realised in the real system.