Modelling and optimisation for sustainable ... - Semantic Scholar
European Journal of Operational Research 164 (2005) 475–490 www.elsevier.com/locate/dsw
O.R. Applications
Modelling and optimisation for sustainable development policy assessment M. Cannon *, B. Kouvaritakis, G. Huang Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, UK Received 9 May 2003; accepted 12 December 2003 Available online 26 February 2004
Abstract Earlier work on sustainable development devised a policy assessment tool that was based on a static optimisation formulation. Key ingredients in sustainable development problems are the presence of random effects and the conflict between different objectives. To accommodate these, the earlier formulation was strongly stochastic and was posed in a multi-objective framework. The purpose of this paper is to consider the extension of the work to a formulation that deploys dynamic optimisation. In particular it is the aim here to use simulations based on a large scale model to derive dynamic rather than static representations, to integrate these into the optimisation scheme and to assess the benefits. 2004 Elsevier B.V. All rights reserved. Keywords: Sustainable development; Modelling; Stochastic systems; Optimisation
1. Introduction Sustainable development is progressively becoming an issue of universal concern and of paramount importance in many aspects of human endeavour. It is therefore not surprising to see an exponential growth in the amount of research effort, yet the unpredictable nature of the problem has meant that most of this research tends to be qualitative and often obtuse. It is of course possible to cast the problem in a probabilistic setting and then perform an optimisation of expected values, but that can often lead to meaningless results given the strongly stochastic nature of the problem. Optimisation is used in this context to obtain the ‘‘best’’ strategy, yet given the many conflicting objectives and political priorities, often the concept of ‘‘best’’ is itself meaningless. Recent work [1,2] proposed a formulation which overcame these problems by considering an optimisation that retained both the stochastic and multi-objective nature of policy making/assessment for sustainable development. The key to this approach was a static model describing the cumulative effects, at the end of a forecast horizon, of a single adjustment in the values of ‘‘actions’’ on ‘‘performance indicators’’.
*
Corresponding author. Tel.: +44-1865-273000/273189; fax: +44-1865-273906. E-mail address: [email protected] (M. Cannon).
0377-2217/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2003.12.018
476
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
Common sense suggests that a dynamic model which anticipates the possibility of future adjustments of the values of actions would introduce extra degrees of freedom to be deployed in the optimisation problem. The viability of such an extension is supported by a feasibility study in [3]. However the concern in [3] was with the salient features of the receding horizon use of dynamic models in general stochastic control problems, which implies a different set of objectives more appropriate to tracking/regulation problems. The aim of this paper is to explore the benefits of dynamic modelling in policy assessment for sustainable development problems relating to Research and Development budget allocation between alternative technologies. The system is defined by a detailed multiplicative economic model, PROMETHEUS [4], comprising hundreds of random coefficients and random exogenous variables; a brief description of PROMETHEUS together with the earlier static optimisation is undertaken in Section 2. As in the feasibility study of [3], it is assumed here that the behaviour of PROMETHEUS can be modelled adequately by linear time invariant Weighting Sequence (WS) models. In order to effect a comparison with the earlier static approach, the appropriate prediction equations are developed in Section 3. This development enables the determination of the extra degrees of freedom that are introduced through the use of dynamics. The prediction equations of Section 3 are deployed in Section 4 to explore the benefits of discounting and to suggest a systematic procedure for the selection of discounting rates. The modelling objective of the paper is undertaken in Section 5 where the earlier assumption concerning the validity of finite length linear timeinvariant WS models is tested by means of Monte Carlo simulations on PROMETHEUS. The end product of these simulations is the computation of the expected values and covariance matrices for the distributions of such WS models. These are used to perform the optimisation of the predicted trajectories for the current and future Research and Development budget allocation and the results thus obtained are compared with those possible through the static optimisation of approach of [1,2]. The paper ends with a set of conclusions in Section 6.
2. Earlier work The system to be considered in this paper is that described by the stochastic multiplicative self-contained energy model PROMETHEUS (for details see [4]), developed by ICCS-NTUA for the European Commission SAPIENT project (Contract no. ENG2-CT1999-0003) and implemented using EVIEWS. It comprises variables that relate to demographic and economic indicators, energy (fuel) consumption, fuel resources and prices, CO2 emissions, technology uptake, and two-factor learning curves describing technology improvement in terms of research and experience gained through application. A summary of the specification for PROMETHEUS [4] is given in Fig. 1 in the form of a flow chart. Exogenous variables, parameters and error terms are modelled as random variables with given distributions derived through extensive use of Monte Carlo simulations and econometric techniques involving a process of variance estimation, normalization, decomposition, and rescaling [4]. Implicit in this is a normality assumption on the error terms which lead to results that contradict economic theory (e.g. parameter non-sign-definiteness) with non-trivial probability (given the large size of variances in PROMETHEUS). The compounded effect of such contradictions could be unacceptable, and for this reason illegal values (together with associated parameter values, whether illegal or not) are rejected, thereby necessitating the performance of further Monte Carlo simulations. PROMETHEUS has been used in the derivation of static models that form the basis of a tool for policy assessment with regard to research and development (R&D) budget allocations. Orthogonal shocks of up to 10% were applied at the beginning of a forecast horizon to the different actions, namely R&D budget allocation in technologies such as Biomass Gasification Gas Turbines (BGGT), Gas Turbine in Combined Cycle (CCGT), Wind Turbines (WIND). Their effects on targets such as cumulative CO2 emissions (CO2 ), or averaged and discounted energy cost (EC), were modelled through equations such as
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
477
PROMETHEUS SPECIFICATION (Summary Flow Chart)
Demographic and Macro Reduced Form Sub-Model
Final Demand Sub-Model
Consumer Prices
Power Generation Sub-Model
Technology Costs
Natural Fossil Fuel Endowments
Primary Prices
Excess Demand Supply Measure
Primary Fuel Demand
Generating Capacity by Technology
2 Factor Learning Curves
Reserves
R&D Expenditure by Technology
TARGETS: Market Impact CO2 Emissions Employment Total En. Cost Security of Supply ……..
Fig. 1. PROMETHEUS specification.
yi ¼ riT x;
ð2:1Þ
where yi is a measure associated with the various targets (e.g. CO2 , EC), ri is a vector of impact coefficients associated with yi and x is the vector of actions (budget allocations), xj , associated with the different technologies. The individual elements of ri , denoted rij , are computed as the ratio of the change in the ith objective measured in specific units to the size of the shock (change in R&D expenditure in euro) applied to the jth action. The above calculation also involves discounting (taken to be 4%) as a means of differentiating between the emphasis on cost/benefit now and the cost/benefit sometime in the future (up to say 30 years ahead). Discounting is common practice in economic and sustainable development modelling and will be applied throughout the sequel. Given the stochastic nature of PROMETHEUS, the impacts ri are random variables whose distribution is taken to be normal: ri Nðqi ; Vi Þ;
ð2:2Þ
where qi , Vi are the mean and covariance matrix of the distribution and are computed on the basis of Monte Carlo simulations (involving approximately 1000 runs each). In reality none of the distributions turn out to be normal, but this assumption will be endorsed because it has a minimal effect on the results of the optimisation problem to be discussed below. Suppose now that a particular target y1 represents cumulative benefit (say at the end of a forecast horizon of 30 years) and the other targets yi , i > 1, represent ()1) times the value of cumulative costs (over the same horizon). Then a sensible policy would aim at the maximization of y1 under the constraint that yi remain above some desired bounds, Ai : max y1 x
subject to yi > Ai ;
ð2:3Þ
478
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
where x here denotes a vector of shocks applied to the actions (budget allocations) at the beginning of the forecast horizon. Clearly all the elements of x must be positive, and the sum of the elements of x must respect a total budget allocation constraint. These two constraints on x are rather obvious and therefore will be assumed throughout the paper (without explicit mention). However, given the random nature of ri , the optimisation in this form could only be applied to the expected values of yj , and this would ignore the significant stochastic character of problems in sustainable development. Instead the problem can be recast in a stochastic framework as [1,2]: max Prfy1 > A1 g x
subject to Prfyi > Ai g P pi
ð2:4Þ
for a set of thresholds Ai and probabilities pi . Note that this formulation preserves the multi-objective nature of the problem in that y1 can be interchanged with any particular yi . Furthermore, exploration of the effect of different choices for Ai , pi provides the flexibility required by a tool for policy assessment. However to remain within a scenario of hedging rather than gambling, all the probabilities have to be chosen to be greater than 0.5; a very desirable consequence of this restriction is that it ensures that the optimisation of (2.4) is convex and therefore has a unique optimum which can be computed reliably and efficiently using Second-Order Conic Programming (SOCP) [8]. The optimisation of (2.4) is in essence ‘‘open-loop’’ in that it only uses information up to current time to improve performance in the future. At subsequent times however, more information will be available and the way that this can be exploited is through a receding horizon application of (2.4). Under this scheme (2.4) is repeated at each new instant of ‘‘real’’ (as opposed to ‘‘prediction’’) time and although the budget allocation is assumed constant over the forecast horizon it will be updated as dictated by the current optimal. This will result in closed-loop performance which will differ from the open-loop optimal, and provides the feedback mechanism required for correction due to model mismatch. Clearly, open-loop optimality does not imply closed-loop optimality, but the latter would require the use of an infinite forecast horizon which, given uncertainty would be meaningless and would result in an infinite-dimensional optimisation problem. The formulation of (2.4) produces useful insights but involves a static optimisation in that budget allocation adjustments are applied only once at the beginning of the forecast horizon and what is being optimised is the cumulative effect of this adjustment at some future points which lies well beyond the transient behaviour of the targets yi . This restriction is apparent from the form of (2.1), which takes no explicit account of the dimension of time, as measured, for example, in terms of number of years within the forecast horizon. It is intuitively obvious that an optimisation of the sort described by (2.4), which considered the same objective function, but which allowed for the possibility that budget allocations can be adjusted in the future (within the forecast horizon), would introduce extra degrees of freedom that could therefore lead to improved optimal solutions. The incorporation of these extra degrees of freedom requires the development of a dynamic model which considers the dependence of on yi ðk þ ijkÞ, the i-steps-ahead prediction for yi at time k, on the current, past and future (predicted) values of x. A preliminary feasibility study [3] speculated that it may be possible to achieve this through the use of Autoregressive Moving Average (ARMA) or Weighting Sequence (WS) models. Of these, the latter is preferable because it is then possible to preserve the assumption of normal distributions for the predicted values of yi . The approach taken in [3] considered the connections that may exist between the receding horizon optimisation strategy as a tool for assessment of sustainable development policy and a general stochastic model predictive control problem. Accordingly there was concern for converting the sustainable development problem into one of set-point tracking, using an equality constraint to guarantee closed-loop stability, and illustrating the efficacy of the proposed scheme in terms of closed-loop responses. The purpose of the current paper is to use PROMETHEUS and test the feasibility of approximating the prediction of suitably discounted dynamics by those of linear time-invariant WS models and to use
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
479
such models in order to explore the benefits of dynamic optimisation. As explained earlier there is no direct relationship between open-loop and closed-loop results, though of course open-loop optimisation should in most cases lead to improvements of closed-loop performance. To avoid this ambiguity, the advantages of dynamic versus static optimisation will be examined entirely in the context of an openloop optimisation. The WS models used in [3] were hypothetical and related to the case of a single action x, in the presence of two yi variables. Here consideration is focused on the assessment of sustainable development policy in connection with two targets: EC and CO2 , and two technologies: CCGT and WIND. 3. Dynamic prediction models and their benefits Given linearity, it has been argued [3] that ARMA and WS models provide convenient representation of system dynamics for the purposes of introducing further degrees of freedom into the optimisation problem of (2.4). Due to the stochastic nature of the PROMETHEUS parameters and exogenous variables, the coefficients of either the ARMA or WS models will be random variables, and this suggests that WS models are to be preferred. The reason for this is that output predictions are by definition random and regression on the predictions therefore involves the multiplication of random predicted values by random model coefficients. By contrast WS models involve the multiplication by random model coefficients of predicted values of inputs that represent deterministic degrees of freedom. The testing of the linearity hypothesis and the computation of WS models for the Sustainable Development Integrated Policy Assessment will be undertaken in Section 5. The aim of this section is to consider the type of prediction equations implied by WS models in the presence of discounting, and to investigate the potential benefits of dynamic modelling versus the static formulation described in Section 2. Thus let an impulsive shock applied to the jth input produce, at successive time instants (measured in years), a change in the ith output response: f gij ð0Þ gij ð1Þ
gij ðqÞ
0
0 g;
where it is implicitly assumed that after q years the response becomes insignificant and can thus be assumed to be zero. As discussed in Section 4, this assumption is made possible through the use of discounting. Then linearity suggests the following dynamic dependence of the output y on the input u: ykþ1 ¼ g0 uk þ g1 uk1 þ þ gq ukq ;
ð3:1Þ
where, for clarity of presentation, we have suppressed the subscripts i; j and instead have used subscripts to indicate time dependencies. This is done with the understanding that similar equations can be written for each i; j pair. However the values generated by PROMETHEUS can be represented by: ykþ1 ¼ g0 uk þ g1 uk1 þ þ gq ukq þ ekþ1 ;
ð3:2Þ
where ekþ1 is introduced to account for model mismatch and for the stochastic nature of PROMETHEUS. Under the assumption that ekþ1 is normally distributed with zero mean and given variance, it is possible to use a least squares procedure to obtain the best estimates for the parameters gi of the WS model. Of course these estimates would themselves be random, normally distributed, and statistical analysis (e.g. [5]) can be deployed to derive the mean and variance: g Nð g; V Þ:
ð3:3Þ
480
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
At time k ¼ 0, the WS model of (3.1) implies a set of prediction equations: 3 3 2 2 3 3 2 2 3 2 g0 u0 þ g1 u1 þ þ gq uq y1 uq u0 u1 7 6 y 7 6 g u þ g u þ þ g u 7 7 6u 6 7 6 0 1 1 0 q qþ1 7 6 2 7 6 6 qþ1 7 6 u1 7 6 u0 7 7 7 6 6 7 7 7 6 6 6 6 y3 7 6 g0 u2 þ g1 u1 þ þ gq uqþ2 7 6 uqþ2 7 6 u2 7 6 u1 7 7 7 6 6 7 7 7 6 6 6 7 6 . 7 6 7 7 6u 6 7 6 .. 7 6 . 7 6 u u 7: 7 7 6 6 6 qþ3 3 2 . 7 ¼ g0 6 7 þ g1 6 6 . 7¼6 7 7 þ þ gq 6 7 7 6 6 7 7 7 6 6 6 u u4 7 u3 7 6 yN 7 6 g0 uN 1 þ g1 uN 2 þ þ gq uN q 7 qþ4 7 6 6 6 7 7 6 6 7 7 7 6 6 6 7 7 6 6 6 uqþ5 7 6 u5 7 6 u4 7 6 yN þ1 7 6 g0 uN þ g1 uN 1 þ þ gq uN qþ1 7 5 5 5 4 4 4 5 5 4 4 .. .. .. .. . . . . . . .
ð3:4Þ
Thus defining wi as the sum of the first i values of y, we have: wN ¼ ðg0 eTN u þ g1 eTN 1 u þ þ gq eTN q u Þ þ ðg1 f1T u þ g2 f2T u þ þ gq fqT u Þ; !
!
!
ð3:5Þ
where u ¼ ½ u0
!
u1
u2
u3
T
uN 1 ;
u ¼ ½ u1
u2
uqþ1
T
ð3:6Þ
and where eTN i is the row vector of 1Õs, with the last i being replaced by 0Õs, and fiT is the vector of 0Õs with the first i being replaced by 1Õs. The above can be written more compactly as wN ¼ gT Uu þ gT Lu ;
ð3:7Þ
!
where U ; L have the form: 2 1 1 1 6 61 1 1 U ¼6. . .. .. 4 .. .. . . 1 1 1
1 .. .
1 .. .
1 0
0 0
.. . .. .
1
3
7 07 ; .. 7 .5 0
2
1 61 L¼6 4 ...
0 1 .. .
.. .
3 0 07 : .. 7 .5
1
1
1
ð3:8Þ
Thus L is a ðq 1Þ ðq 1Þ lower triangular matrix of 1Õs, whereas U comprises the columns of L taken in reverse order, preceded by a q ðN qÞ matrix of 1Õs. Before exploring the available degrees of freedom, it is necessary to ensure that y reaches a steady state by the end of the forecast horizon, otherwise the maximization of wN would provide an unreliable indicator of ‘‘benefit/cost’’. This can only be achieved by fixing all predicted inputs after the first m ðm 6 N qÞ to a constant value, which will be denoted by uss , the subscript indicating steady state. Under such circumstances, the prediction equation for wN becomes: 2 3 1 N m 61 N m 17 6 7 T 7; wN ¼ g ðMx þ Lu Þ; M ¼ 6 .. 6 .. 7 4. 5 . ð3:9Þ 1 N mq m1 X x ¼ ½ x1 x2 T ; x1 ¼ ui ; x2 ¼ uss : i¼0
Thus re-introducing the subscripts i; j, it is possible to write separate prediction equation for each of the cumulative outputs as:
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
wi ðk þ N þ 1jkÞ ¼
m X
gjT ðMxj þ Lu ðkÞÞ;
xj ¼ ½ x1j
T
x2j ;
x1j ¼
j
j¼1
m1 X
uj ðk þ ljkÞ;
481
x2j ¼ uj ðssjkÞ;
l¼0
ð3:10Þ where the denotation ‘‘k þ ljk’’ has been used to emphasise the fact that the values above refer to predictions based on information up to time k; for simplicity only it has been assumed that the values for q; m are the same for each i; j pair. The above can be written more compactly as: wi ðk þ N þ 1jkÞ ¼ hTi x þ bi ðkÞ; hTi ¼ gi1T M
T gim M ;
x ¼ xT1
xTm
T
;
bi ðkÞ ¼
q X j¼1
gijT Lu ðkÞ:
ð3:11Þ
j
This form of the prediction equations is similar to that of (2.1) with two important differences noted in the remarks below. Remark 3.1. Eq. (2.1) is linear in the degrees of freedom, whereas, on account of dynamics which involve past values of inputs (accounted for by bi ðkÞ), Eq. (3.11) is affine rather than linear; this difference does not apply to the case of ‘‘open-loop optimisation’’ considered in this paper, because what is being considered is incremental behaviour and this implies zero initial conditions (i.e. bi ðkÞ ¼ 0). Remark 3.2. The vector of degrees of freedom in (3.11) is twice as long as that of (2.4). This observation may appear surprising given the fact that each input is now allowed to assume different values over the whole of the forecast horizon; this would imply (3.11) contains significantly more than just twice the number of degrees of freedom present in (2.4). However, as explained above the requirement to reach steady state within the forecast horizon removes q degrees of freedom and replaces them by just one, ui ðssjkÞ, for each of the inputs. In addition, given the cumulative nature of wi , their values are affected by the sum (rather than the individual choice) of the first m predicted values of each input. As indicated above, the vector of WS parameters for each pair i; j is random, and the least squares identification procedure implies a normal distribution as indicated in (3.3). As a consequence, the vectors hi will also be random and normally distributed and simple analysis can be used to compute the mean and covariance matrix: hi ; V i Þ ð3:12Þ hi Nð from the mean and covariance matrices for the vector of parameters for each i; j pair, or more directly (as is done in Section 5) from Monte Carlo simulations performed on PROMETHEUS directly. Therefore it is possible to formulate the ‘‘dynamic’’ equivalent of the ‘‘static’’ optimisation problem of (2.4) as: max Prfw1 ðk þ N þ 1jkÞ > A1 g subject to Prfwi ðk þ N þ 1jkÞ > Ai g P pi x
ð3:13Þ
or more conveniently: A1 hT1 x min pffiffiffiffiffiffiffiffiffiffiffi x xT V1 x
Ai hTi x ffi 6 ai ; subject to pffiffiffiffiffiffiffiffiffiffi xT Vi x
ð3:14Þ
where ai is such that Prfz 6 ai g ¼ pi , for z Nð0; 1Þ. In Section 5 suitable models are derived (together with their mean and covariance matrices) describing the dynamic effects of two inputs, u1 , u2 (constituting actions on the R&D allocation between two technologies, CCGT and WIND) on two outputs, y1 , y2 (measuring performance of two indicators, one
482
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
concerning CO2 and the other EC). First however, consideration will be given to the effect of discounting with particular reference to the justification of using Finite Impulse Response (FIR) models, namely WS models with a finite number of non-zero parameters. For clarity, once again the i; j subscripts will be dropped and instead subscripts will denote time-dependence; the understanding being that the whole procedure can be repeated for each i; j pair, using the appropriate WS model fgij ðkÞg.
4. Discounting and FIRs Sustainability refers to the strategy of aiming for the attainment of development at current time while retaining some potential for development for future generations, hence the selection of a 30 year forecast horizon. However the present and the future are not given equal weighting, it being deemed that benefit/ cost now is far more important than benefit/cost 30 years in the future. The mechanism for shifting the emphasis from the future to the current time is discounting, according to which the i-step-ahead prediction for y is scaled by qi , where q is a positive number strictly less than one (often in the range of 0.9–0.96). Thus denoting the discounted output values by ~y , the predictions of (3.4) become: 32 3 2 3 2 ~y1 g0 0 0 0 u0 7 6 ~y2 7 6 qg1 6 qg0 0 0 7 7 6 u1 7 6 7 6 2 2 7 6 ~y3 7 6 q2 g2 6 q g1 q g0 0 7 6 u2 7 ð4:1Þ 7; 6 7¼6 6 .. 7 6 .. .. 76 .. 7 .. .. 4 . 5 4 . . 54 . 5 . . uN 1 qN 1 gN 1 qN 1 gN 2 qN 1 gN 3 qN 1 g0 ~yN where, as explained earlier, due to the open-loop nature of the result considered in this paper the effects of the past have been set equal to zero. The prediction equations above are not in convenient form in that they cannot be generated by a linear time-invariant WS model. A convenient solution to this, which concurrently makes sense within the Sustainable Development framework, is to discount the inputs as well as the outputs. In this case the above becomes: 32 3 2 3 2 ~y1 g0 0 0 0 ~u0 7 6 ~y2 7 6 qg1 6 g0 0 0 7 76 ~u1 7 6 7 6 76 ~u2 7 6 ~y3 7 6 q2 g2 qg g 0 1 0 ð4:2Þ 76 7 6 7¼6 6 .. 7 6 .. .. 76 .. 7 .. .. 5 5 4 . 5 4 4 . . . . . ~yN
qN 1 gN 1
qN 2 gN 2
qN 3 gN 3
g0
~uN 1
thereby yielding the prediction equations for a linear time-invariant system with WS model: g~0 g~1 g~2 ¼ g0 qg1 q2 g2 :
ð4:3Þ
This forms the basis for the justification that the PROMETHEUS dynamics can be represented by WS models with a finite number, q, of non-zero terms: the discounting sequence f1; q; q2 ; . . .g converges to zero, and q can always be chosen to be small enough so that the discounted WS above will itself converge to zero. Under such circumstances, the discounted WS can be truncated (incurring a sufficiently small approximation error) after a finite number of terms, i.e. q, where the choice of q clearly depends on the choice of q. Thus consideration of the requirement that steady state should be reached within the forecast horizon, N , implies that q 6 N m 6 N 1 and this in turn implies an upper bound on q. Of course, if there is no concern about whether or not the steady state is reached within the forecast horizon, then q could be taken to be as long as the forecast horizon ðq ¼ N Þ, and in this case the role of discounting is not so vital in the definition of FIRÕs. This is particularly relevant to CO2 emissions, which are normally not discounted. Mixing discounted inputs with undiscounted outputs would contradict the use of a time-invariant WS
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
483
prediction model; instead one would have to resort to multi-model prediction equations [7]. For simplicity, in this paper all inputs and outputs will be discounted.
5. Computation of dynamic models and optimisation It was assumed in the earlier sections that the system behaviour (i.e. the behaviour of PROMETHEUS) can be represented (to within reasonable accuracy) by linear dynamic models. This hypothesis was tested for a variety of inputs that remained within a total budget constraint corresponding to a 10% increase, and was found to be true. It is not claimed here that there was no evidence of non-linear behaviour and that one could not do better through the use of non-linear models (e.g. Hammerstein models), but rather that for actions within the given budget allocation, the system was sufficiently weakly non-linear so as to justify the use of linear WS models. Fig. 2 provides the results of a typical test where an input sequence varying randomly from 0 to 1 over the entire forecast horizon, scaled so as to be on the limit of the total budget allocation, was applied to CCGT; for convenience, such input sequences will be referred to as FPRBS (Feasible Pseudo Random Binary Sequences). The same input sequence after scaling by a factor 0.75 was also applied to CCGT, and the ratios of the corresponding responses (CO2 and EC) computed at each instant of forecast time were plotted, and the whole procedure was repeated 5 times (in total) with respect to both CCGT and WIND inputs. From Fig. 2, it is apparent that these ratios were close to the input scaling of 0.75. The spikes in the fourth plot are outliers corresponding to zero crossings when the ratios involve
y12
y11
0.8
0.76
0.79
0.74
0.78 0.72 0.77 0.7 0.76 0.68
0.75 0.74
0
10
20
30
0.66
0
10
y21
20
30
20
30
y22 4
0.8
3 2
0.78
1 0.76 0 0.74
0
10
20
30
-1
0
Fig. 2. Linearity test.
10
484
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
division by values close to zero. The same effect (but to a lesser extent) is observed for the first few forecast instants and is due to immediate transients during which some of the responses are close to zero. The WS model in Eqs. (4.2) and (4.3) can be rewritten in the form: 3 3 2 3 2 2 ~y1 g~0 ~ u0 0 0 0 6 ~y2 7 6 ~ 6 g~1 7 ~ u0 0 0 7 7 7 6 7 6 u1 6 7 7 6 ~y3 7 6 ~ ~ u2 0 7; ~h ¼ 6 u1 u0 h þ ~e; D ¼ 6 ~ ð5:1Þ 6 7 ¼ D~ 6 g~2 7; 6 .. 7 6 .. 6 .. 7 .. .. 7 .. 4 . 5 4 . 4 . 5 . . 5 . ~ uN 1 ~ uN 2 ~ uN 3 ~u0 ~yN g~N 1 where ~e is a vector of residuals introduced to account for the mismatch between the behaviour of the system (PROMETHEUS) and the model. This equation could be used as the basis of a least-squares identification procedure. However, given that PROMETHEUS provides the system behaviour up to 30 instants of time (this being the forecast horizon), it was necessary to repeat the experiment over a large number of times (in fact 100 times) using FPRBS input sequences, and then to concatenate the results so as to produce a single identification problem with a number of observations (30 · 100) which (as required) far exceeds the number of parameters to be identified ðq 6 30Þ. Normally the optimal truncation level q would be obtained through the application of procedures such as the Akaike criterion [6], but here it was found that the optimum was in excess of the permissible upper bound of q < 30. Fig. 3 shows the optimal results for q ¼ 20 (circles) and for q ¼ 30 (dashed line); from these plots it is apparent that truncation after q ¼ 20 is perfectly feasible. A comparison undertaken in Fig. 4 of the WS models for q ¼ 20 (solid lines) and the PROMETHEUS
0.5
g11
x 10–3
1
0
g12
x 10–3
0
–0.5 –1 –1 –2 –1.5 –3
–2 –2.5
2
0
10
20
30
–4
g21
x 10–4
1.5
0
10
20
30
20
30
g22
x 10–4
1 0
0.5 0
–2 –0.5 –1
–4 0
10
20
30
–1.5
0
10
Fig. 3. Identified WS models for q ¼ 20 (circles) and q ¼ 30 (dashed line).
5
1
0
0
–5
–1
y12(k)
y11(k)
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
–10
–2
–15
–3
–20
–4
–25
0
10
20
–5
30
485
0
10
k
20
30
20
30
k
1
0.2
0
0.1
y22(k)
y21(k)
–1 –2
0 –0.1
–3 –0.2
–4 –5
0
10
20 k
30
–0.3
0
10 k
Fig. 4. Impulse responses: WS model (solid line), PROMETHEUS (dashed line).
impulse responses (dashed lines) obtained upon the application of orthogonal shocks, shows that with the exception of the pair i ¼ 1, j ¼ 2, agreement is good, the differences reflecting the non-linear and timevarying behaviour of PROMETHEUS. A more meaningful comparison would be to contrast the response of the WS models (dashed lines) against that of PROMETHEUS (solid lines) for the same input sequences and this is undertaken in Fig. 5 where two randomly selected FPRBS are applied, one for CCGT and the other for WIND. Agreement here is seen to be very good and the plots provide convincing evidence as to the suitability of the WS approximations. The formulation of the optimisation problem of (3.14) requires the expected values and covariance matrices for hi . The former can be computed from (3.11) when vectors giT are replaced by the vectors of the elements of the identified WS shown in Fig. 3. Alternatively, both the expected values and covariance matrices for hi can be computed from Monte Carlo simulations of (3.11) on PROMETHEUS; as explained earlier, bi ðkÞ is to be taken to be zero. It is noted that such simulations will model the incremental effect of the predicted shocks ui ðkj0Þ and as such will include the effects of the random PROMETHEUS coefficients, but not those of the exogenous variables. This is fortuitous because, depending on the size of the total budget, the effects of exogenous variables could dominate the behaviour of the various performance indicators thereby obscuring any beneficial effects that the policy of optimising the predicted allocations ui ðkj0Þ might have. It is emphasised that the purpose of the tool being developed is policy assessment, not the minimization of the effect random exogenous variables (which after all could only be effected through the use of feedback and would therefore concern a receding horizon implementation of the optimisation strategy). The results of the Monte Carlo simulations are summarised below. 9:7249e 3 2:9901e 5 2:3382e 4 T M g~11 N ; ; ð5:2Þ 8:0158e 2 2:3382e 4 1:8340e 3
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490 0
0
–50
–20
y12(k)
y11(k)
486
–100
–150
–200
–40
–60
0
10
20
–80
30
0
10
k 0
20
30
0.5
y22(k)
y21(k)
30
1
–5
–10
–15
–20
20 k
0 –0.5 –1
0
10
20
30
–1.5
0
10
k
k
Fig. 5. Responses to FPRBS: WS model (solid line), PROMETHEUS (dashed line).
2:6021e 2 1:9301e 4 M g~12 N ; 1:9682e 1 1:4977e 3
1:4977e 3 1:1875e 2
1:1327e 3 2:6544e 7 M g~21 N ; 9:4616e 3 2:1270e 6
2:1270e 6 1:7212e 5
1:2517e 4 1:3508e 6 M g~22 N ; 6:6997e 4 9:6789e 6
9:6789e 6 7:0022e 5
T
T
T
;
ð5:3Þ
;
ð5:4Þ
:
ð5:5Þ
Introducing (5.2)–(5.4) into the static optimisation of (2.4), which can be implemented through the formulation of (3.14) provided that that x1i and x2i are constrained to be equal to one another for all i, results in the optimal discounted input sequences shown in Fig. 6. The data for this figure, and for Figs. 7–9 was: A1 ¼ 2000, p1 ¼ 0:7, and A2 ¼ 200 with p2 being the probability that is being maximised. The corresponding output sequences are shown in Fig. 7, where the results are shown as obtained from the WS models (solid lines) as well as through PROMETHEUS (dashed lines) when this is used in its deterministic (expected values) form. The respective results for the dynamic formulation of (3.14) are shown in Figs. 8 and 9 and can be seen to be vastly different, despite the fact that the improvement in the probability that is being maximised is somewhat modest (the static optimal was p2 ¼ 0:5646, and the dynamic optimal p2 ¼ 0:5825, representing an improvement of 3.2%). The significant difference in input sequences is also reflected in the resulting output sequences, and indeed the expected cumulative effect on EC, as computed by PROMETHEUS, is )208.3 for the static formulation and )233.7 for the dynamic; these figures represent reductions (with respect to a baseline) on energy costs, and it is clear that the dynamic results outperform (by about 12%) the static results.
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
487
12000 10000
u1(k)
8000 6000 4000 2000 0
0
5
10
15
20
25
30
0
5
10
15
20
25
30
2500 2000
u2(k)
1500 1000 500 0
k Fig. 6. Inputs for static optimisation.
0
y1(k)
–50 –100 –150 –200
0
5
10
15
20
25
30
0
5
10
15
20
25
30
0 –2
y2(k)
–4 –6 –8 –10 –12 –14
k Fig. 7. Outputs for static optimisation: WS (solid), PROMETHEUS (dashed).
The modest improvement in maximised probability can be made a great deal more significant (rising by about 12% from 0.5646 for the static to 0.6311 for the dynamic) through the use of all degrees of freedom that are available in the dynamic formulation. However this improvement is achieved through input sequences that become very active towards the end of the forecast horizon. From a sustainable development
488
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
3.5
x 104
3
u1(k)
2.5 2 1.5 1 0.5 0
0
5
10
15
20
25
30
0
5
10
15
20
25
30
10000
u2(k)
8000 6000 4000 2000 0
k Fig. 8. Inputs for dynamic optimisation.
50 0
y1(k)
–50 –100 –150 –200 –250 –300
0
5
10
15
20
25
30
0
5
10
15
20
25
30
5 0
y2(k)
–5 –10 –15 –20 –25
k Fig. 9. Outputs for dynamic optimisation: WS (solid), PROMETHEUS (dashed).
perspective, such behaviour may be less useful and for this reason the corresponding input and output sequences will not be included in this paper. Fig. 10 summarises the benefits possible in terms of the maximised probability p2 as the threshold A2 of the cumulative sum of y2 is allowed to vary; the results for the static, dynamic with 2 degrees of freedom and dynamic with all the degrees of freedom are plotted in dashed, solid and dotted lines respectively. Clearly
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
489
1 0.95 0.9 0.85
max p2
0.8 0.75 0.7 0.65 0.6 0.55 0.5 –250
–200
–150
–100
–50
0
50
100
A2 Fig. 10. The effect of A2 on the optimal values of p2 .
the effect of A2 is significant and this is in direct contrast to the effects of p1 and A1 which are small; for this reason no further plots have been included. Whereas the variation of the maximised p2 probability shown in Fig. 10 is continuous, the effect of this variation on the optimising input trajectories for the dynamic formulation with 2 degrees of freedom is
u1(av,k
O.R. Applications
Modelling and optimisation for sustainable development policy assessment M. Cannon *, B. Kouvaritakis, G. Huang Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, UK Received 9 May 2003; accepted 12 December 2003 Available online 26 February 2004
Abstract Earlier work on sustainable development devised a policy assessment tool that was based on a static optimisation formulation. Key ingredients in sustainable development problems are the presence of random effects and the conflict between different objectives. To accommodate these, the earlier formulation was strongly stochastic and was posed in a multi-objective framework. The purpose of this paper is to consider the extension of the work to a formulation that deploys dynamic optimisation. In particular it is the aim here to use simulations based on a large scale model to derive dynamic rather than static representations, to integrate these into the optimisation scheme and to assess the benefits. 2004 Elsevier B.V. All rights reserved. Keywords: Sustainable development; Modelling; Stochastic systems; Optimisation
1. Introduction Sustainable development is progressively becoming an issue of universal concern and of paramount importance in many aspects of human endeavour. It is therefore not surprising to see an exponential growth in the amount of research effort, yet the unpredictable nature of the problem has meant that most of this research tends to be qualitative and often obtuse. It is of course possible to cast the problem in a probabilistic setting and then perform an optimisation of expected values, but that can often lead to meaningless results given the strongly stochastic nature of the problem. Optimisation is used in this context to obtain the ‘‘best’’ strategy, yet given the many conflicting objectives and political priorities, often the concept of ‘‘best’’ is itself meaningless. Recent work [1,2] proposed a formulation which overcame these problems by considering an optimisation that retained both the stochastic and multi-objective nature of policy making/assessment for sustainable development. The key to this approach was a static model describing the cumulative effects, at the end of a forecast horizon, of a single adjustment in the values of ‘‘actions’’ on ‘‘performance indicators’’.
*
Corresponding author. Tel.: +44-1865-273000/273189; fax: +44-1865-273906. E-mail address: [email protected] (M. Cannon).
0377-2217/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2003.12.018
476
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
Common sense suggests that a dynamic model which anticipates the possibility of future adjustments of the values of actions would introduce extra degrees of freedom to be deployed in the optimisation problem. The viability of such an extension is supported by a feasibility study in [3]. However the concern in [3] was with the salient features of the receding horizon use of dynamic models in general stochastic control problems, which implies a different set of objectives more appropriate to tracking/regulation problems. The aim of this paper is to explore the benefits of dynamic modelling in policy assessment for sustainable development problems relating to Research and Development budget allocation between alternative technologies. The system is defined by a detailed multiplicative economic model, PROMETHEUS [4], comprising hundreds of random coefficients and random exogenous variables; a brief description of PROMETHEUS together with the earlier static optimisation is undertaken in Section 2. As in the feasibility study of [3], it is assumed here that the behaviour of PROMETHEUS can be modelled adequately by linear time invariant Weighting Sequence (WS) models. In order to effect a comparison with the earlier static approach, the appropriate prediction equations are developed in Section 3. This development enables the determination of the extra degrees of freedom that are introduced through the use of dynamics. The prediction equations of Section 3 are deployed in Section 4 to explore the benefits of discounting and to suggest a systematic procedure for the selection of discounting rates. The modelling objective of the paper is undertaken in Section 5 where the earlier assumption concerning the validity of finite length linear timeinvariant WS models is tested by means of Monte Carlo simulations on PROMETHEUS. The end product of these simulations is the computation of the expected values and covariance matrices for the distributions of such WS models. These are used to perform the optimisation of the predicted trajectories for the current and future Research and Development budget allocation and the results thus obtained are compared with those possible through the static optimisation of approach of [1,2]. The paper ends with a set of conclusions in Section 6.
2. Earlier work The system to be considered in this paper is that described by the stochastic multiplicative self-contained energy model PROMETHEUS (for details see [4]), developed by ICCS-NTUA for the European Commission SAPIENT project (Contract no. ENG2-CT1999-0003) and implemented using EVIEWS. It comprises variables that relate to demographic and economic indicators, energy (fuel) consumption, fuel resources and prices, CO2 emissions, technology uptake, and two-factor learning curves describing technology improvement in terms of research and experience gained through application. A summary of the specification for PROMETHEUS [4] is given in Fig. 1 in the form of a flow chart. Exogenous variables, parameters and error terms are modelled as random variables with given distributions derived through extensive use of Monte Carlo simulations and econometric techniques involving a process of variance estimation, normalization, decomposition, and rescaling [4]. Implicit in this is a normality assumption on the error terms which lead to results that contradict economic theory (e.g. parameter non-sign-definiteness) with non-trivial probability (given the large size of variances in PROMETHEUS). The compounded effect of such contradictions could be unacceptable, and for this reason illegal values (together with associated parameter values, whether illegal or not) are rejected, thereby necessitating the performance of further Monte Carlo simulations. PROMETHEUS has been used in the derivation of static models that form the basis of a tool for policy assessment with regard to research and development (R&D) budget allocations. Orthogonal shocks of up to 10% were applied at the beginning of a forecast horizon to the different actions, namely R&D budget allocation in technologies such as Biomass Gasification Gas Turbines (BGGT), Gas Turbine in Combined Cycle (CCGT), Wind Turbines (WIND). Their effects on targets such as cumulative CO2 emissions (CO2 ), or averaged and discounted energy cost (EC), were modelled through equations such as
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
477
PROMETHEUS SPECIFICATION (Summary Flow Chart)
Demographic and Macro Reduced Form Sub-Model
Final Demand Sub-Model
Consumer Prices
Power Generation Sub-Model
Technology Costs
Natural Fossil Fuel Endowments
Primary Prices
Excess Demand Supply Measure
Primary Fuel Demand
Generating Capacity by Technology
2 Factor Learning Curves
Reserves
R&D Expenditure by Technology
TARGETS: Market Impact CO2 Emissions Employment Total En. Cost Security of Supply ……..
Fig. 1. PROMETHEUS specification.
yi ¼ riT x;
ð2:1Þ
where yi is a measure associated with the various targets (e.g. CO2 , EC), ri is a vector of impact coefficients associated with yi and x is the vector of actions (budget allocations), xj , associated with the different technologies. The individual elements of ri , denoted rij , are computed as the ratio of the change in the ith objective measured in specific units to the size of the shock (change in R&D expenditure in euro) applied to the jth action. The above calculation also involves discounting (taken to be 4%) as a means of differentiating between the emphasis on cost/benefit now and the cost/benefit sometime in the future (up to say 30 years ahead). Discounting is common practice in economic and sustainable development modelling and will be applied throughout the sequel. Given the stochastic nature of PROMETHEUS, the impacts ri are random variables whose distribution is taken to be normal: ri Nðqi ; Vi Þ;
ð2:2Þ
where qi , Vi are the mean and covariance matrix of the distribution and are computed on the basis of Monte Carlo simulations (involving approximately 1000 runs each). In reality none of the distributions turn out to be normal, but this assumption will be endorsed because it has a minimal effect on the results of the optimisation problem to be discussed below. Suppose now that a particular target y1 represents cumulative benefit (say at the end of a forecast horizon of 30 years) and the other targets yi , i > 1, represent ()1) times the value of cumulative costs (over the same horizon). Then a sensible policy would aim at the maximization of y1 under the constraint that yi remain above some desired bounds, Ai : max y1 x
subject to yi > Ai ;
ð2:3Þ
478
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
where x here denotes a vector of shocks applied to the actions (budget allocations) at the beginning of the forecast horizon. Clearly all the elements of x must be positive, and the sum of the elements of x must respect a total budget allocation constraint. These two constraints on x are rather obvious and therefore will be assumed throughout the paper (without explicit mention). However, given the random nature of ri , the optimisation in this form could only be applied to the expected values of yj , and this would ignore the significant stochastic character of problems in sustainable development. Instead the problem can be recast in a stochastic framework as [1,2]: max Prfy1 > A1 g x
subject to Prfyi > Ai g P pi
ð2:4Þ
for a set of thresholds Ai and probabilities pi . Note that this formulation preserves the multi-objective nature of the problem in that y1 can be interchanged with any particular yi . Furthermore, exploration of the effect of different choices for Ai , pi provides the flexibility required by a tool for policy assessment. However to remain within a scenario of hedging rather than gambling, all the probabilities have to be chosen to be greater than 0.5; a very desirable consequence of this restriction is that it ensures that the optimisation of (2.4) is convex and therefore has a unique optimum which can be computed reliably and efficiently using Second-Order Conic Programming (SOCP) [8]. The optimisation of (2.4) is in essence ‘‘open-loop’’ in that it only uses information up to current time to improve performance in the future. At subsequent times however, more information will be available and the way that this can be exploited is through a receding horizon application of (2.4). Under this scheme (2.4) is repeated at each new instant of ‘‘real’’ (as opposed to ‘‘prediction’’) time and although the budget allocation is assumed constant over the forecast horizon it will be updated as dictated by the current optimal. This will result in closed-loop performance which will differ from the open-loop optimal, and provides the feedback mechanism required for correction due to model mismatch. Clearly, open-loop optimality does not imply closed-loop optimality, but the latter would require the use of an infinite forecast horizon which, given uncertainty would be meaningless and would result in an infinite-dimensional optimisation problem. The formulation of (2.4) produces useful insights but involves a static optimisation in that budget allocation adjustments are applied only once at the beginning of the forecast horizon and what is being optimised is the cumulative effect of this adjustment at some future points which lies well beyond the transient behaviour of the targets yi . This restriction is apparent from the form of (2.1), which takes no explicit account of the dimension of time, as measured, for example, in terms of number of years within the forecast horizon. It is intuitively obvious that an optimisation of the sort described by (2.4), which considered the same objective function, but which allowed for the possibility that budget allocations can be adjusted in the future (within the forecast horizon), would introduce extra degrees of freedom that could therefore lead to improved optimal solutions. The incorporation of these extra degrees of freedom requires the development of a dynamic model which considers the dependence of on yi ðk þ ijkÞ, the i-steps-ahead prediction for yi at time k, on the current, past and future (predicted) values of x. A preliminary feasibility study [3] speculated that it may be possible to achieve this through the use of Autoregressive Moving Average (ARMA) or Weighting Sequence (WS) models. Of these, the latter is preferable because it is then possible to preserve the assumption of normal distributions for the predicted values of yi . The approach taken in [3] considered the connections that may exist between the receding horizon optimisation strategy as a tool for assessment of sustainable development policy and a general stochastic model predictive control problem. Accordingly there was concern for converting the sustainable development problem into one of set-point tracking, using an equality constraint to guarantee closed-loop stability, and illustrating the efficacy of the proposed scheme in terms of closed-loop responses. The purpose of the current paper is to use PROMETHEUS and test the feasibility of approximating the prediction of suitably discounted dynamics by those of linear time-invariant WS models and to use
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
479
such models in order to explore the benefits of dynamic optimisation. As explained earlier there is no direct relationship between open-loop and closed-loop results, though of course open-loop optimisation should in most cases lead to improvements of closed-loop performance. To avoid this ambiguity, the advantages of dynamic versus static optimisation will be examined entirely in the context of an openloop optimisation. The WS models used in [3] were hypothetical and related to the case of a single action x, in the presence of two yi variables. Here consideration is focused on the assessment of sustainable development policy in connection with two targets: EC and CO2 , and two technologies: CCGT and WIND. 3. Dynamic prediction models and their benefits Given linearity, it has been argued [3] that ARMA and WS models provide convenient representation of system dynamics for the purposes of introducing further degrees of freedom into the optimisation problem of (2.4). Due to the stochastic nature of the PROMETHEUS parameters and exogenous variables, the coefficients of either the ARMA or WS models will be random variables, and this suggests that WS models are to be preferred. The reason for this is that output predictions are by definition random and regression on the predictions therefore involves the multiplication of random predicted values by random model coefficients. By contrast WS models involve the multiplication by random model coefficients of predicted values of inputs that represent deterministic degrees of freedom. The testing of the linearity hypothesis and the computation of WS models for the Sustainable Development Integrated Policy Assessment will be undertaken in Section 5. The aim of this section is to consider the type of prediction equations implied by WS models in the presence of discounting, and to investigate the potential benefits of dynamic modelling versus the static formulation described in Section 2. Thus let an impulsive shock applied to the jth input produce, at successive time instants (measured in years), a change in the ith output response: f gij ð0Þ gij ð1Þ
gij ðqÞ
0
0 g;
where it is implicitly assumed that after q years the response becomes insignificant and can thus be assumed to be zero. As discussed in Section 4, this assumption is made possible through the use of discounting. Then linearity suggests the following dynamic dependence of the output y on the input u: ykþ1 ¼ g0 uk þ g1 uk1 þ þ gq ukq ;
ð3:1Þ
where, for clarity of presentation, we have suppressed the subscripts i; j and instead have used subscripts to indicate time dependencies. This is done with the understanding that similar equations can be written for each i; j pair. However the values generated by PROMETHEUS can be represented by: ykþ1 ¼ g0 uk þ g1 uk1 þ þ gq ukq þ ekþ1 ;
ð3:2Þ
where ekþ1 is introduced to account for model mismatch and for the stochastic nature of PROMETHEUS. Under the assumption that ekþ1 is normally distributed with zero mean and given variance, it is possible to use a least squares procedure to obtain the best estimates for the parameters gi of the WS model. Of course these estimates would themselves be random, normally distributed, and statistical analysis (e.g. [5]) can be deployed to derive the mean and variance: g Nð g; V Þ:
ð3:3Þ
480
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
At time k ¼ 0, the WS model of (3.1) implies a set of prediction equations: 3 3 2 2 3 3 2 2 3 2 g0 u0 þ g1 u1 þ þ gq uq y1 uq u0 u1 7 6 y 7 6 g u þ g u þ þ g u 7 7 6u 6 7 6 0 1 1 0 q qþ1 7 6 2 7 6 6 qþ1 7 6 u1 7 6 u0 7 7 7 6 6 7 7 7 6 6 6 6 y3 7 6 g0 u2 þ g1 u1 þ þ gq uqþ2 7 6 uqþ2 7 6 u2 7 6 u1 7 7 7 6 6 7 7 7 6 6 6 7 6 . 7 6 7 7 6u 6 7 6 .. 7 6 . 7 6 u u 7: 7 7 6 6 6 qþ3 3 2 . 7 ¼ g0 6 7 þ g1 6 6 . 7¼6 7 7 þ þ gq 6 7 7 6 6 7 7 7 6 6 6 u u4 7 u3 7 6 yN 7 6 g0 uN 1 þ g1 uN 2 þ þ gq uN q 7 qþ4 7 6 6 6 7 7 6 6 7 7 7 6 6 6 7 7 6 6 6 uqþ5 7 6 u5 7 6 u4 7 6 yN þ1 7 6 g0 uN þ g1 uN 1 þ þ gq uN qþ1 7 5 5 5 4 4 4 5 5 4 4 .. .. .. .. . . . . . . .
ð3:4Þ
Thus defining wi as the sum of the first i values of y, we have: wN ¼ ðg0 eTN u þ g1 eTN 1 u þ þ gq eTN q u Þ þ ðg1 f1T u þ g2 f2T u þ þ gq fqT u Þ; !
!
!
ð3:5Þ
where u ¼ ½ u0
!
u1
u2
u3
T
uN 1 ;
u ¼ ½ u1
u2
uqþ1
T
ð3:6Þ
and where eTN i is the row vector of 1Õs, with the last i being replaced by 0Õs, and fiT is the vector of 0Õs with the first i being replaced by 1Õs. The above can be written more compactly as wN ¼ gT Uu þ gT Lu ;
ð3:7Þ
!
where U ; L have the form: 2 1 1 1 6 61 1 1 U ¼6. . .. .. 4 .. .. . . 1 1 1
1 .. .
1 .. .
1 0
0 0
.. . .. .
1
3
7 07 ; .. 7 .5 0
2
1 61 L¼6 4 ...
0 1 .. .
.. .
3 0 07 : .. 7 .5
1
1
1
ð3:8Þ
Thus L is a ðq 1Þ ðq 1Þ lower triangular matrix of 1Õs, whereas U comprises the columns of L taken in reverse order, preceded by a q ðN qÞ matrix of 1Õs. Before exploring the available degrees of freedom, it is necessary to ensure that y reaches a steady state by the end of the forecast horizon, otherwise the maximization of wN would provide an unreliable indicator of ‘‘benefit/cost’’. This can only be achieved by fixing all predicted inputs after the first m ðm 6 N qÞ to a constant value, which will be denoted by uss , the subscript indicating steady state. Under such circumstances, the prediction equation for wN becomes: 2 3 1 N m 61 N m 17 6 7 T 7; wN ¼ g ðMx þ Lu Þ; M ¼ 6 .. 6 .. 7 4. 5 . ð3:9Þ 1 N mq m1 X x ¼ ½ x1 x2 T ; x1 ¼ ui ; x2 ¼ uss : i¼0
Thus re-introducing the subscripts i; j, it is possible to write separate prediction equation for each of the cumulative outputs as:
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
wi ðk þ N þ 1jkÞ ¼
m X
gjT ðMxj þ Lu ðkÞÞ;
xj ¼ ½ x1j
T
x2j ;
x1j ¼
j
j¼1
m1 X
uj ðk þ ljkÞ;
481
x2j ¼ uj ðssjkÞ;
l¼0
ð3:10Þ where the denotation ‘‘k þ ljk’’ has been used to emphasise the fact that the values above refer to predictions based on information up to time k; for simplicity only it has been assumed that the values for q; m are the same for each i; j pair. The above can be written more compactly as: wi ðk þ N þ 1jkÞ ¼ hTi x þ bi ðkÞ; hTi ¼ gi1T M
T gim M ;
x ¼ xT1
xTm
T
;
bi ðkÞ ¼
q X j¼1
gijT Lu ðkÞ:
ð3:11Þ
j
This form of the prediction equations is similar to that of (2.1) with two important differences noted in the remarks below. Remark 3.1. Eq. (2.1) is linear in the degrees of freedom, whereas, on account of dynamics which involve past values of inputs (accounted for by bi ðkÞ), Eq. (3.11) is affine rather than linear; this difference does not apply to the case of ‘‘open-loop optimisation’’ considered in this paper, because what is being considered is incremental behaviour and this implies zero initial conditions (i.e. bi ðkÞ ¼ 0). Remark 3.2. The vector of degrees of freedom in (3.11) is twice as long as that of (2.4). This observation may appear surprising given the fact that each input is now allowed to assume different values over the whole of the forecast horizon; this would imply (3.11) contains significantly more than just twice the number of degrees of freedom present in (2.4). However, as explained above the requirement to reach steady state within the forecast horizon removes q degrees of freedom and replaces them by just one, ui ðssjkÞ, for each of the inputs. In addition, given the cumulative nature of wi , their values are affected by the sum (rather than the individual choice) of the first m predicted values of each input. As indicated above, the vector of WS parameters for each pair i; j is random, and the least squares identification procedure implies a normal distribution as indicated in (3.3). As a consequence, the vectors hi will also be random and normally distributed and simple analysis can be used to compute the mean and covariance matrix: hi ; V i Þ ð3:12Þ hi Nð from the mean and covariance matrices for the vector of parameters for each i; j pair, or more directly (as is done in Section 5) from Monte Carlo simulations performed on PROMETHEUS directly. Therefore it is possible to formulate the ‘‘dynamic’’ equivalent of the ‘‘static’’ optimisation problem of (2.4) as: max Prfw1 ðk þ N þ 1jkÞ > A1 g subject to Prfwi ðk þ N þ 1jkÞ > Ai g P pi x
ð3:13Þ
or more conveniently: A1 hT1 x min pffiffiffiffiffiffiffiffiffiffiffi x xT V1 x
Ai hTi x ffi 6 ai ; subject to pffiffiffiffiffiffiffiffiffiffi xT Vi x
ð3:14Þ
where ai is such that Prfz 6 ai g ¼ pi , for z Nð0; 1Þ. In Section 5 suitable models are derived (together with their mean and covariance matrices) describing the dynamic effects of two inputs, u1 , u2 (constituting actions on the R&D allocation between two technologies, CCGT and WIND) on two outputs, y1 , y2 (measuring performance of two indicators, one
482
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
concerning CO2 and the other EC). First however, consideration will be given to the effect of discounting with particular reference to the justification of using Finite Impulse Response (FIR) models, namely WS models with a finite number of non-zero parameters. For clarity, once again the i; j subscripts will be dropped and instead subscripts will denote time-dependence; the understanding being that the whole procedure can be repeated for each i; j pair, using the appropriate WS model fgij ðkÞg.
4. Discounting and FIRs Sustainability refers to the strategy of aiming for the attainment of development at current time while retaining some potential for development for future generations, hence the selection of a 30 year forecast horizon. However the present and the future are not given equal weighting, it being deemed that benefit/ cost now is far more important than benefit/cost 30 years in the future. The mechanism for shifting the emphasis from the future to the current time is discounting, according to which the i-step-ahead prediction for y is scaled by qi , where q is a positive number strictly less than one (often in the range of 0.9–0.96). Thus denoting the discounted output values by ~y , the predictions of (3.4) become: 32 3 2 3 2 ~y1 g0 0 0 0 u0 7 6 ~y2 7 6 qg1 6 qg0 0 0 7 7 6 u1 7 6 7 6 2 2 7 6 ~y3 7 6 q2 g2 6 q g1 q g0 0 7 6 u2 7 ð4:1Þ 7; 6 7¼6 6 .. 7 6 .. .. 76 .. 7 .. .. 4 . 5 4 . . 54 . 5 . . uN 1 qN 1 gN 1 qN 1 gN 2 qN 1 gN 3 qN 1 g0 ~yN where, as explained earlier, due to the open-loop nature of the result considered in this paper the effects of the past have been set equal to zero. The prediction equations above are not in convenient form in that they cannot be generated by a linear time-invariant WS model. A convenient solution to this, which concurrently makes sense within the Sustainable Development framework, is to discount the inputs as well as the outputs. In this case the above becomes: 32 3 2 3 2 ~y1 g0 0 0 0 ~u0 7 6 ~y2 7 6 qg1 6 g0 0 0 7 76 ~u1 7 6 7 6 76 ~u2 7 6 ~y3 7 6 q2 g2 qg g 0 1 0 ð4:2Þ 76 7 6 7¼6 6 .. 7 6 .. .. 76 .. 7 .. .. 5 5 4 . 5 4 4 . . . . . ~yN
qN 1 gN 1
qN 2 gN 2
qN 3 gN 3
g0
~uN 1
thereby yielding the prediction equations for a linear time-invariant system with WS model: g~0 g~1 g~2 ¼ g0 qg1 q2 g2 :
ð4:3Þ
This forms the basis for the justification that the PROMETHEUS dynamics can be represented by WS models with a finite number, q, of non-zero terms: the discounting sequence f1; q; q2 ; . . .g converges to zero, and q can always be chosen to be small enough so that the discounted WS above will itself converge to zero. Under such circumstances, the discounted WS can be truncated (incurring a sufficiently small approximation error) after a finite number of terms, i.e. q, where the choice of q clearly depends on the choice of q. Thus consideration of the requirement that steady state should be reached within the forecast horizon, N , implies that q 6 N m 6 N 1 and this in turn implies an upper bound on q. Of course, if there is no concern about whether or not the steady state is reached within the forecast horizon, then q could be taken to be as long as the forecast horizon ðq ¼ N Þ, and in this case the role of discounting is not so vital in the definition of FIRÕs. This is particularly relevant to CO2 emissions, which are normally not discounted. Mixing discounted inputs with undiscounted outputs would contradict the use of a time-invariant WS
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
483
prediction model; instead one would have to resort to multi-model prediction equations [7]. For simplicity, in this paper all inputs and outputs will be discounted.
5. Computation of dynamic models and optimisation It was assumed in the earlier sections that the system behaviour (i.e. the behaviour of PROMETHEUS) can be represented (to within reasonable accuracy) by linear dynamic models. This hypothesis was tested for a variety of inputs that remained within a total budget constraint corresponding to a 10% increase, and was found to be true. It is not claimed here that there was no evidence of non-linear behaviour and that one could not do better through the use of non-linear models (e.g. Hammerstein models), but rather that for actions within the given budget allocation, the system was sufficiently weakly non-linear so as to justify the use of linear WS models. Fig. 2 provides the results of a typical test where an input sequence varying randomly from 0 to 1 over the entire forecast horizon, scaled so as to be on the limit of the total budget allocation, was applied to CCGT; for convenience, such input sequences will be referred to as FPRBS (Feasible Pseudo Random Binary Sequences). The same input sequence after scaling by a factor 0.75 was also applied to CCGT, and the ratios of the corresponding responses (CO2 and EC) computed at each instant of forecast time were plotted, and the whole procedure was repeated 5 times (in total) with respect to both CCGT and WIND inputs. From Fig. 2, it is apparent that these ratios were close to the input scaling of 0.75. The spikes in the fourth plot are outliers corresponding to zero crossings when the ratios involve
y12
y11
0.8
0.76
0.79
0.74
0.78 0.72 0.77 0.7 0.76 0.68
0.75 0.74
0
10
20
30
0.66
0
10
y21
20
30
20
30
y22 4
0.8
3 2
0.78
1 0.76 0 0.74
0
10
20
30
-1
0
Fig. 2. Linearity test.
10
484
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
division by values close to zero. The same effect (but to a lesser extent) is observed for the first few forecast instants and is due to immediate transients during which some of the responses are close to zero. The WS model in Eqs. (4.2) and (4.3) can be rewritten in the form: 3 3 2 3 2 2 ~y1 g~0 ~ u0 0 0 0 6 ~y2 7 6 ~ 6 g~1 7 ~ u0 0 0 7 7 7 6 7 6 u1 6 7 7 6 ~y3 7 6 ~ ~ u2 0 7; ~h ¼ 6 u1 u0 h þ ~e; D ¼ 6 ~ ð5:1Þ 6 7 ¼ D~ 6 g~2 7; 6 .. 7 6 .. 6 .. 7 .. .. 7 .. 4 . 5 4 . 4 . 5 . . 5 . ~ uN 1 ~ uN 2 ~ uN 3 ~u0 ~yN g~N 1 where ~e is a vector of residuals introduced to account for the mismatch between the behaviour of the system (PROMETHEUS) and the model. This equation could be used as the basis of a least-squares identification procedure. However, given that PROMETHEUS provides the system behaviour up to 30 instants of time (this being the forecast horizon), it was necessary to repeat the experiment over a large number of times (in fact 100 times) using FPRBS input sequences, and then to concatenate the results so as to produce a single identification problem with a number of observations (30 · 100) which (as required) far exceeds the number of parameters to be identified ðq 6 30Þ. Normally the optimal truncation level q would be obtained through the application of procedures such as the Akaike criterion [6], but here it was found that the optimum was in excess of the permissible upper bound of q < 30. Fig. 3 shows the optimal results for q ¼ 20 (circles) and for q ¼ 30 (dashed line); from these plots it is apparent that truncation after q ¼ 20 is perfectly feasible. A comparison undertaken in Fig. 4 of the WS models for q ¼ 20 (solid lines) and the PROMETHEUS
0.5
g11
x 10–3
1
0
g12
x 10–3
0
–0.5 –1 –1 –2 –1.5 –3
–2 –2.5
2
0
10
20
30
–4
g21
x 10–4
1.5
0
10
20
30
20
30
g22
x 10–4
1 0
0.5 0
–2 –0.5 –1
–4 0
10
20
30
–1.5
0
10
Fig. 3. Identified WS models for q ¼ 20 (circles) and q ¼ 30 (dashed line).
5
1
0
0
–5
–1
y12(k)
y11(k)
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
–10
–2
–15
–3
–20
–4
–25
0
10
20
–5
30
485
0
10
k
20
30
20
30
k
1
0.2
0
0.1
y22(k)
y21(k)
–1 –2
0 –0.1
–3 –0.2
–4 –5
0
10
20 k
30
–0.3
0
10 k
Fig. 4. Impulse responses: WS model (solid line), PROMETHEUS (dashed line).
impulse responses (dashed lines) obtained upon the application of orthogonal shocks, shows that with the exception of the pair i ¼ 1, j ¼ 2, agreement is good, the differences reflecting the non-linear and timevarying behaviour of PROMETHEUS. A more meaningful comparison would be to contrast the response of the WS models (dashed lines) against that of PROMETHEUS (solid lines) for the same input sequences and this is undertaken in Fig. 5 where two randomly selected FPRBS are applied, one for CCGT and the other for WIND. Agreement here is seen to be very good and the plots provide convincing evidence as to the suitability of the WS approximations. The formulation of the optimisation problem of (3.14) requires the expected values and covariance matrices for hi . The former can be computed from (3.11) when vectors giT are replaced by the vectors of the elements of the identified WS shown in Fig. 3. Alternatively, both the expected values and covariance matrices for hi can be computed from Monte Carlo simulations of (3.11) on PROMETHEUS; as explained earlier, bi ðkÞ is to be taken to be zero. It is noted that such simulations will model the incremental effect of the predicted shocks ui ðkj0Þ and as such will include the effects of the random PROMETHEUS coefficients, but not those of the exogenous variables. This is fortuitous because, depending on the size of the total budget, the effects of exogenous variables could dominate the behaviour of the various performance indicators thereby obscuring any beneficial effects that the policy of optimising the predicted allocations ui ðkj0Þ might have. It is emphasised that the purpose of the tool being developed is policy assessment, not the minimization of the effect random exogenous variables (which after all could only be effected through the use of feedback and would therefore concern a receding horizon implementation of the optimisation strategy). The results of the Monte Carlo simulations are summarised below. 9:7249e 3 2:9901e 5 2:3382e 4 T M g~11 N ; ; ð5:2Þ 8:0158e 2 2:3382e 4 1:8340e 3
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490 0
0
–50
–20
y12(k)
y11(k)
486
–100
–150
–200
–40
–60
0
10
20
–80
30
0
10
k 0
20
30
0.5
y22(k)
y21(k)
30
1
–5
–10
–15
–20
20 k
0 –0.5 –1
0
10
20
30
–1.5
0
10
k
k
Fig. 5. Responses to FPRBS: WS model (solid line), PROMETHEUS (dashed line).
2:6021e 2 1:9301e 4 M g~12 N ; 1:9682e 1 1:4977e 3
1:4977e 3 1:1875e 2
1:1327e 3 2:6544e 7 M g~21 N ; 9:4616e 3 2:1270e 6
2:1270e 6 1:7212e 5
1:2517e 4 1:3508e 6 M g~22 N ; 6:6997e 4 9:6789e 6
9:6789e 6 7:0022e 5
T
T
T
;
ð5:3Þ
;
ð5:4Þ
:
ð5:5Þ
Introducing (5.2)–(5.4) into the static optimisation of (2.4), which can be implemented through the formulation of (3.14) provided that that x1i and x2i are constrained to be equal to one another for all i, results in the optimal discounted input sequences shown in Fig. 6. The data for this figure, and for Figs. 7–9 was: A1 ¼ 2000, p1 ¼ 0:7, and A2 ¼ 200 with p2 being the probability that is being maximised. The corresponding output sequences are shown in Fig. 7, where the results are shown as obtained from the WS models (solid lines) as well as through PROMETHEUS (dashed lines) when this is used in its deterministic (expected values) form. The respective results for the dynamic formulation of (3.14) are shown in Figs. 8 and 9 and can be seen to be vastly different, despite the fact that the improvement in the probability that is being maximised is somewhat modest (the static optimal was p2 ¼ 0:5646, and the dynamic optimal p2 ¼ 0:5825, representing an improvement of 3.2%). The significant difference in input sequences is also reflected in the resulting output sequences, and indeed the expected cumulative effect on EC, as computed by PROMETHEUS, is )208.3 for the static formulation and )233.7 for the dynamic; these figures represent reductions (with respect to a baseline) on energy costs, and it is clear that the dynamic results outperform (by about 12%) the static results.
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
487
12000 10000
u1(k)
8000 6000 4000 2000 0
0
5
10
15
20
25
30
0
5
10
15
20
25
30
2500 2000
u2(k)
1500 1000 500 0
k Fig. 6. Inputs for static optimisation.
0
y1(k)
–50 –100 –150 –200
0
5
10
15
20
25
30
0
5
10
15
20
25
30
0 –2
y2(k)
–4 –6 –8 –10 –12 –14
k Fig. 7. Outputs for static optimisation: WS (solid), PROMETHEUS (dashed).
The modest improvement in maximised probability can be made a great deal more significant (rising by about 12% from 0.5646 for the static to 0.6311 for the dynamic) through the use of all degrees of freedom that are available in the dynamic formulation. However this improvement is achieved through input sequences that become very active towards the end of the forecast horizon. From a sustainable development
488
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
3.5
x 104
3
u1(k)
2.5 2 1.5 1 0.5 0
0
5
10
15
20
25
30
0
5
10
15
20
25
30
10000
u2(k)
8000 6000 4000 2000 0
k Fig. 8. Inputs for dynamic optimisation.
50 0
y1(k)
–50 –100 –150 –200 –250 –300
0
5
10
15
20
25
30
0
5
10
15
20
25
30
5 0
y2(k)
–5 –10 –15 –20 –25
k Fig. 9. Outputs for dynamic optimisation: WS (solid), PROMETHEUS (dashed).
perspective, such behaviour may be less useful and for this reason the corresponding input and output sequences will not be included in this paper. Fig. 10 summarises the benefits possible in terms of the maximised probability p2 as the threshold A2 of the cumulative sum of y2 is allowed to vary; the results for the static, dynamic with 2 degrees of freedom and dynamic with all the degrees of freedom are plotted in dashed, solid and dotted lines respectively. Clearly
M. Cannon et al. / European Journal of Operational Research 164 (2005) 475–490
489
1 0.95 0.9 0.85
max p2
0.8 0.75 0.7 0.65 0.6 0.55 0.5 –250
–200
–150
–100
–50
0
50
100
A2 Fig. 10. The effect of A2 on the optimal values of p2 .
the effect of A2 is significant and this is in direct contrast to the effects of p1 and A1 which are small; for this reason no further plots have been included. Whereas the variation of the maximised p2 probability shown in Fig. 10 is continuous, the effect of this variation on the optimising input trajectories for the dynamic formulation with 2 degrees of freedom is
u1(av,k
