Sensitivity and Covariance in Stochastic Complementarity Problems with an Application to Natural Gas Markets
SSensitivity and Covariance in StochasticComplementarity Problems with an Application toNorth American Natural Gas Markets
Sriram Sankaranarayanan a, ∗ , Felipe Feijoo b , Sauleh Siddiqui a,c,d a Department of Civil Engineering, Johns Hopkins University b Pontificia Universidad Cat´olica de Valpara´ıso, Chile c Department of Applied Mathematics and Statistics, Johns Hopkins University d German Institute for Economic Research (DIW Berlin)
Abstract
We provide an efficient method to approximate the covariance between decisionvariables and uncertain parameters in solutions to a general class of stochasticnonlinear complementarity problems. We also develop a sensitivity metric toquantify uncertainty propagation by determining the change in the varianceof the output due to a change in the variance of an input parameter. Thecovariance matrix of the solution variables quantifies the uncertainty in theoutput and pairs correlated variables and parameters. The sensitivity metrichelps in identifying the parameters that cause maximum fluctuations in theoutput. The method developed in this paper optimizes the use of gradientsand matrix multiplications which makes it particularly useful for large-scaleproblems. Having developed this method, we extend the deterministic versionof the North American Natural Gas Model (NANGAM), to incorporate effectsdue to uncertainty in the parameters of the demand function, supply function,infrastructure costs, and investment costs. We then use the sensitivity metricsto identify the parameters that impact the equilibrium the most.
Keywords:
Stochastic Programming, Large scale optimization,Complementarity problems, Approximation methods ∗ Corresponding Author: [email protected]
Preprint submitted to Elsevier a r X i v : . [ m a t h . O C ] N ov . Introduction Complementarity models arise naturally out of various real life problems.A rigorous survey of their application is available in [19]. Authors in [1, 11,18, 28, 37, 40] use complementarity problems to model markets from a gametheoretic perspective [4, 44], where the complementarity conditions typicallyarise between the marginal profit and the quantity produced by the producer.In the field of mechanics, they typically arise in the context of frictional contactproblems [34], where there is a complementarity relation between the frictionalforce between a pair of surfaces and the distance of separation between them.With a wide range of applications, understanding the characteristics of solutionsto complementarity problems becomes important for advancing the field. In thispaper, we focus on studying the characteristics of solutions to complementarityproblems under uncertainty.The behavior of a solution to a complementarity problem with random pa-rameters was first addressed in [25], where such problems were referred to asstochastic complementarity problems (SCP). Authors in [9, 14, 23, 30, 42] de-fine various formulations of SCP for different applications and have devisedalgorithms to solve the problem. Authors in [35] compute confidence intervalsfor solution of the expected value formulation of the problem, however theydo not have efficient methods to find the second-order statistics for large-scalecomplementarity problems.Large-scale problems, those with over 10,000 decision variables and uncertainparameters arise naturally out of detailed market models and there is consider-able interest in studying, understanding and solving such models. For example,[10] discuss a case of urban drainage system with large number of variables. [45]discuss a case of deciding under large-scale nuclear emergencies. In line withthe area of application used in this paper, [22] discuss a case of an energy modelwith large number of variables and parameters. Naturally, developing methodsto solve such large-scale problems gained interest. Authors in [32, 36] discussvarious tools ranging from mathematical techniques (decomposition based) to2omputational techniques (parallel processing) for solving large-scale optimiza-tion problems. [39] uses an approximate algorithm for a large-scale Markov de-cision process to optimize production and distribution systems. In this paper,we do not present a new method to solve stochastic complementarity problems,but an efficient algorithm to generate second-order information that is flexibleenough to be coupled with any existing algorithm that provides a first-ordersolution.The objective of this paper is to efficiently obtain second-order statistical in-formation about solution vectors of large-scale stochastic complementarity prob-lems. This gives us information about variability of the equilibrium obtainedby solving a nonlinear complementarity problem (NCP) and the correlation be-tween various variables in the solution. Authors in [29] and [6] provide examplesin the area of clinical pathways and ecology respectively, about the utility ofunderstanding the variance of the solution in addition to the mean. They alsoshow that a knowledge of variance aids better understanding and planning ofthe system. [3] emphasize the necessity to understand covariance as a wholerather than individual variances by quantifying “the loss incurred on ignoringcorrelations” in a stochastic programming model.In addition, we also introduce a sensitivity metric which quantifies the changein uncertainty in the output due to a perturbation in the variance of uncertaininput parameters. This helps us to directly compare input parameters by theamount of uncertainty they propagate to the solution.In attaining the above objectives, the most computationally expensive stepis to solve a system of linear equations. We choose approximation methodsover analytical methods, integration, or Monte Carlo simulation because of thecomputational hurdle involved while implementing those methods for large-scaleproblems. The method we describe in this paper achieves the following: • The most expensive step has to be performed just once, irrespective of thecovariance of the input parameters. Once the linear system of equationsis solved, for each given covariance scenario, we only perform two matrix3ultiplications. • Approximating the covariance matrix and getting a sensitivity metric canbe obtained by solving the above mentioned linear system just once.The methods developed in this paper can also be used for nonlinear opti-mization problems with linear equality constraints. We prove stronger resultson error bounds for special cases of quadratic programming.Having developed this method, we apply it to a large-scale stochastic naturalgas model for North America, an extension of the deterministic model developedin [18] and determine the covariance of the solution variables. We then proceedto identify the parameters which have the greatest impact on the solution. APython class for efficiently storing and operating on sparse arrays of dimensiongreater than two is created. This is useful for working with high-dimensionalproblems which have an inherent sparse structure in the gradients.We divide the paper as follows. Section 2 formulates the problem and men-tions the assumptions used in the paper. It then develops the algorithm usedto approximate the solution covariance and provides proofs for bounding theerror. Section 3 develops a framework to quantify the sensitivity of the solutionto each of the random variables. Section 4 shows how the result can be appliedfor certain optimization problems with equality constraints. Having obtainedthe theoretical results, section 5 gives an example of a oligopoly where thismethod can be applied and compares the computational time of the approxima-tion method with a Monte-Carlo method showing the performance improvementfor large-scale problems. Section 6 describes the Natural Gas Model to whichthe said method is applied. Section 7 discusses the possible enhancements forthe model and its limitations in the current form.
2. Approximation of covariance
For the rest of the paper, all bold quantities are vectors. A subscript i forthose quantities refer to the i -th component of the vector in Cartesian represen-tation. 4 .1. Definitions We define a complementarity problem and a stochastic complementarityproblem which are central to the results obtained in this paper. We use ageneral definition of complementarity problems and stochastic complementarityproblems as stated below.
Definition 1. [17] Given F : R n × m (cid:55)→ R n , and parameters θ ∈ R m , the parametrized nonlinear complementarity problem (NCP) is to find x ∈ R n suchthat K (cid:51) x ⊥ F ( x ; θ ) ∈ K ∗ (2.1)where K ∗ , the dual cone of K is defined as K ∗ = (cid:8) x ∈ R n : v T x ≥ ∀ v ∈ K (cid:9) (2.2) Definition 2.
Given a cone K ∈ R n a random function F : K × Ω (cid:55)→ R n , thestochastic complementarity problem (SCP) is to find x ∈ R n such that K (cid:51) x ⊥ E F ( x ; ω ) ∈ K ∗ (2.3)We assume that we can explicitly evaluate the expectation in (2.3) using itsfunctional form, and that the SCP can be solved using an existing algorithm.We now make assumptions on the form of K in (2.1). This form of K helps in establishing an equivalence between a complementarity problem and aminimization problem which is key to derive the approximation method in thispaper. Assumption 1. K in (2.1) is a Cartesian product of half spaces and full spaces, i.e., for some I ⊆ { , , . . . , n } K = { x ∈ R n : x i ≥ i ∈ I} (2.4)We now propose a lemma about the form of the dual cone of K to understandthe special form that it has. This will help us convert the complementarityproblem into an unconstrained minimization problem. Lemma 3.
The dual cone K ∗ of the set assumed in assumption 1 is K ∗ = K (cid:48) = (cid:26) x ∈ R n : x i ≥ if i ∈ I x i = 0 if i (cid:54)∈ I (cid:27) (2.5) Proof.
Check Appendix A (cid:4) .2. Preliminaries for approximation In this subsection, we prove two preliminary results. Firstly, we prove ourability to pose an NCP as an unconstrained minimization problem. Then weprove results on twice continuous differentiability of the objective function, thusenabling us to use the rich literature available for smooth unconstrained min-imization. Following that, propositions 5 and 6 help in achieving the formerwhile proposition 7 along with its corollaries help us in achieving the latter.We now define C-functions, which are central to pose the complementarityproblem into an unconstrained optimization problem. The equivalent formu-lation as an unconstrained optimization problem assists us in developing thealgorithm.
Definition 4. [17, Pg. 72] A function ψ : R (cid:55)→ R is a C-function when ψ ( x, y ) = 0 ⇔ x ≥ y ≥ xy = 0 (2.6)We consider the following commonly used C-functions. ψ F B ( x, y ) = (cid:112) x + y − x − y (2.7) ψ min ( x, y ) = min( x, y ) (2.8)Under our assumptions on K , the following two propositions establish theequivalence of the complementarity problem and an unconstrained minimizationproblem. Proposition 5.
Suppose assumption 1 holds. Then every solution x ∗ ( θ ) of theparameterized complementarity problem in (2.1) , is a global minimum of thefollowing function f ( x ; θ ) , Φ i ( x , θ ; F ) = (cid:26) F i ( x , θ ) if i (cid:54)∈ I ψ i ( x i , F i ( x , θ )) if i ∈ I (2.9) f ( x ; θ ) = 12 (cid:107) Φ( x ; θ ; F ) (cid:107) (2.10) with an objective value 0, for some set of not necessarily identical C-functions ψ i .Proof. Check Appendix A (cid:4)
Proposition 6.
Suppose assumption 1 holds. If a solution to the problem in (2.1) exists and x ∗ ( θ ) is an unconstrained global minimizer of f ( x ; θ ) defined in (2.10) , then x ∗ ( θ ) solves the complementarity problem in (2.1) . roof. Check Appendix A (cid:4)
Now given a function F , and a set K which satisfies assumption 1, and asolution of the NCP x ∗ (ˆ θ ) for some fixed θ = ˆ θ , we define a vector valuedfunction Φ : R n × m (cid:55)→ R n component-wise as follows.Φ i ( x , θ ; F ) = F i ( x , θ ) if i (cid:54)∈ I ψ ( x i , F i ( x , θ )) if i ∈ Z ψ ( x i , F i ( x , θ )) otherwise (2.11) f ( x ; θ ) = 12 (cid:107) Φ( x ; θ ; F ) (cid:107) (2.12) Z = (cid:110) i ∈ I : x ∗ i (ˆ θ ) = F i ( x ∗ (ˆ θ ); ˆ θ ) = 0 (cid:111) (2.13)Note that if ψ is a C-function, ψ is also a C-function since ψ = 0 ⇐⇒ ψ = 0.We observe from propositions 5 and 6 that minimizing f ( x ; θ ) over x is equivalentto solving the NCP in (2.1).Now we assume conditions on the smoothness of F so that the solution to aperturbed problem is sufficiently close to the original solution. Assumption 2. F ( x ; θ ) is twice continuously differentiable in x and θ over anopen set containing K .Given that the rest of the analysis is on the function f ( x ; θ ) defined in (2.12),we prove that a sufficiently smooth F and a suitable ψ ensure a sufficientlysmooth f . Proposition 7.
With assumption 2 holding, we state f ( x ; ˆ θ ) defined as in (2.12) is a twice continuously differentiable at x satisfying f ( x ; ˆ θ ) = 0 for any C func-tion ψ provided ψ is twice differentiable at (cid:8) ( a, b ) ∈ R : ψ ( a, b ) = 0 (cid:9) \{ (0 , } with a finitederivative and finite second derivative. ψ vanishes sufficiently fast near the origin. i.e., lim ( a,b ) → (0 , ψ ( a, b ) ∂ ψ ( a, b ) ∂a∂b = 0 (2.14) Proof.
Given that f is a sum of squares, it is sufficient to prove each termindividually is twice continuously differentiable to prove the theorem. Alsosince we are only interested where f vanishes, it is sufficient to prove the aboveproperty for each term where it vanishes.7onsider terms from i (cid:54)∈ I . Since F i is twice continuously differentiable, F i istwice continuously differentiable too.Consider the case i ∈ Z . This means i ∈ I and x ∗ i (ˆ θ ) = F i (cid:16) x ∗ i (ˆ θ ) , ˆ θ (cid:17) = 0.These contribute a ψ term to f . With the notation ψ i ≡ ψ ( x i , F i ( x ; θ )), and δ ij = 1 ⇐⇒ i = j and 0 otherwise we clearly have, ∂ψ i ∂ x j = 4 ψ i (cid:18) ∂ψ i ∂a δ ij + ∂ψ i ∂b ∂ F i ∂ x j (cid:19) = 0 (2.15) ∂ ψ i ∂ x j x k = 12 ψ i (cid:18) ∂ψ i ∂a δ ij + ∂ψ i ∂b ∂ F i ∂ x j (cid:19) (cid:18) ∂ψ i ∂a δ ik + ∂ψ i ∂b ∂ F i ∂ x k (cid:19) + 4 ψ i (cid:18) ∂ ψ i ∂a δ ij δ ik + ∂ ψ i ∂a∂b ∂ F i ∂ x k + other terms (cid:19) (2.16)= 0 (2.17)For the third case, i ∈ I \ Z , we have ∂ψ i ∂ x j = 2 ψ i (cid:18) ∂ψ i ∂a δ ij + ∂ψ i ∂b ∂ F i ∂ x j (cid:19) = 0 (2.18) ∂ ψ i ∂ x j x k = (cid:18) ∂ψ i ∂a δ ij + ∂ψ i ∂b ∂ F i ∂ x j (cid:19) + 2 ψ i (cid:18) ∂ ψ i ∂a δ ij δ ik + ∂ ψ i ∂a∂b ∂ F i ∂ x k + other terms (cid:19) (2.19)= ∂ψ i ∂a δ ij + ∂ψ i ∂b ∂ F i ∂ x j (2.20)Continuity of f at the points of interest follow the continuity of the individualterms at the points. (cid:4) The following corollaries show the existence of C-functions ψ which satisfythe hypothesis of proposition 7. Corollary 8.
With assumption 2 holding, for the choice of C-function ψ = ψ min defined in (2.8) , the function f is twice continuously differentiable at itszeros. Corollary 9.
With assumption 2 holding, for the choice of C-function ψ = ψ F B defined in (2.7) , the function f is twice continuously differentiable. We now define an isolated solution to a problem and assume that the problemof interest has this property. This is required to ensure that our approximationis well defined.
Definition 10. [38] A minimum x ∗ of a problem is said to be an isolated min-imum , if there is a neighborhood B ( x ∗ ; (cid:15) ) of x ∗ , where x ∗ is the only minimumof the problem. 8 a) An example of a function where theglobal minimum x = 0 is a non-isolatedsolution (b) The intuition behind our approxi-mation for finding where ∇ f ( x , θ ) = 0under a small perturbation Figure 1: Intuitions
A counter-example for an isolated minimum is shown on Fig. 1a. It is a plotof the function f ( x ) = 5 x + x sin (cid:18) x (cid:19) (2.21)and the global minimum at x = 0 is not an isolated minimum as we can confirmthat any open interval around x = 0 has other minimum contained in it. Unlikethis case, in this paper, we assume that if we obtain a global-minimum of f , thenit is an isolated minimum. The existence of a neighborhood B ( x ∗ ; (cid:15) ) as requiredin definition 10 protects the approximation method from returning to a localminimum in the neighborhood for sufficiently small perturbations in problemparameters. Assumption 3.
For some fixed value of θ = ˆ θ , there exists a known solution x ∗ (ˆ θ ) such that it is an isolated global minimum of f ( x ; θ ). This subsection achieves two primary results as follows. We now proposealgorithm 1 to approximate the covariance of the output given a covariancematrix for the input parameters. Following that Theorem 11 gives a mathemat-ical proof that the algorithm indeed approximates the covariance matrix. Thesecond key result is Theorem 12 which bounds the error in the approximation.9iven a deterministic shift in a parameters value, the sensitivity of the solu-tion has been studied in [7, 20] for mathematical programming problems (e.g.,nonlinear or large-scale programs). Authors in [8] used a perturbation tech-nique to study the sensitivities in calculus of variations. We aim to build onthe research in [7, 20] by proposing an approximation approach for large-scalestochastic complementarity problems. In particular, the sensitivity analysis inChapter 3 in [20] motivates the method developed in this paper, after convertingthe complementarity problem into an unconstrained optimization problem withsufficient smoothness properties. Chapters 4 and 5 in [20] provide good contextfor computational approaches related to the ideas of Chapter 3. Our approachis amenable to extensions towards approximations for large-scale models, as westudy the sensitivity of the variance of the solution, to a perturbation in the vari-ance of input random parameters. We invite readers to refer to these works onfurther context in standard optimization problems. Chapter 2 of [20] providesinteresting examples on how small perturbations can lead to large changes insolutions. This is good context in situations where Assumptions 2 (and Corol-laries 8 and 9) and 3 in our paper do not hold. The research in [7] providesgood context for standard optimization problems, including situations where werelax assumptions on smoothness and active constraints.The intuition behind the approximation is shown on Fig. 1b and can besummarized as follows. Having posed the NCP as an unconstrained minimiza-tion of a function f , we now approximate the change in the solution due toa perturbation of the parameters. Keeping the smoothness properties of f inmind, we say that the gradient of f vanishes at the solution before any pertur-bation. Following the random perturbation of parameters, we approximate thenew value of the gradient of f at the old solution. Then we compute the stepto be taken so that the gradient at the new point vanishes. We formalize thisidea in Theorem 11 and use that to build Algorithm 1 with some features forincreased efficiency. The analysis builds on Dini’s implicit function theorem [33]for deterministic perturbations and extends the results to predict covarianceunder random perturbations. 10 lgorithm 1 Approximating CovarianceSolve the complementarity problem in (2.1) for the mean value of θ = ˆ θ , or solvethe stochastic complementarity problem in (2.3) and calibrate the value of theparameters θ = ˆ θ for this solution. Call this solution as x ∗ . Choose a tolerancelevel τ . Evaluate F ∗ ← F ( x ∗ ; ˆ θ ), G ij ← ∂ F i ( x ∗ ;ˆ θ ) ∂ x j , L i,j ← ∂ F i ( x ∗ ;ˆ θ ) ∂θ j . Choose a C-function ψ such that the conditions in proposition 7 are satisfied. Define the function ψ a ( a, b ) = ∂ψ ( a,b ) ∂a , ψ b ( a, b ) = ∂ψ ( a,b ) ∂b . Find the set of indices Z = { z ∈ I : | x ∗ z | = | F ∗ z | ≤ τ } . Define M ij ← G ij if i (cid:54)∈ I i ∈ Z ψ a ( x ∗ i , F ∗ i ) δ ij + ψ b ( x ∗ i , F ∗ i ) G ij otherwise (2.22)where δ ij = 1 if i = j and 0 otherwise. Define N ij ← L ij if i (cid:54)∈ I i ∈ Z ψ b ( x ∗ i , F ∗ i ) L ij otherwise (2.23) Solve the linear systems of equations for T . MT = N (2.24)If M is non singular, we have a unique solution. If not, a least squaresolution or a solution obtained by calculating the Moore Penrose Pseudoinverse [26] can be used. Given C , a covariance matrix of the input random parameters, θ ( ω ), return C ∗ ← T CT T . Theorem 11.
Algorithm 1 generates Taylor’s first-order approximation for thechange in solution for a perturbation in parameters and computes the covarianceof the solution for a complementarity problem with uncertain parameters withsmall variances.Proof.
Consider the function f ( x ; θ ). From Theorem 5, x ∗ ≡ x ∗ (ˆ θ ) minimizesthis function for θ = ˆ θ . From proposition 7, we have f ( x ; θ ) is twice continuouslydifferentiable at all its zeros. Thus we have, ∇ x f ( x ∗ ; ˆ θ ) = 0 (2.25)Now suppose the parameters ˆ θ are perturbed by ∆ θ , then the above gradientcan be written using the mean value theorem and then approximated up to the11rst order as follows. ∇ x f ( x ∗ (ˆ θ ) , ˆ θ + ∆ θ ) = ∇ x f ( x ∗ (ˆ θ ) , ˆ θ )+ ∇ ˆ θ ∇ x f ( x ∗ (ˆ θ ) , (cid:101) θ )∆ θ (2.26) ∇ x f ( x ∗ (ˆ θ ) , ˆ θ + ∆ θ ) − ∇ x f ( x ∗ (ˆ θ ) , ˆ θ ) ≈ J ∆ θ (2.27)where, (cid:101) θ ∈ [ θ, θ + ∆ θ ] (2.28) J ij = [ ∇ ˆ θ ∇ x f ( x ∗ (ˆ θ ) , ˆ θ )] ij (2.29)= ∂ [ ∇ x f ( x ∗ ; ˆ θ )] i ∂ ˆ θ j (2.30)Since J ∆ θ is not guaranteed to be 0, we might have to alter x to bring thegradient back to zero. i.e., we need ∆ x such that ∇ x f ( x ∗ (ˆ θ ) + ∆ x , ˆ θ + ∆ θ ) = .But by the mean value theorem, ∇ x f ( x ∗ (ˆ θ ) + ∆ x , ˆ θ + ∆ θ ) = ∇ x f ( x ∗ (ˆ θ ) , ˆ θ + ∆ θ )+ ∇ x f ( (cid:101) x , ˆ θ + ∆ θ )∆ x (2.31) ≈ J ∆ θ + ∇ x f ( (cid:101) x , ˆ θ )∆ x (2.32) ≈ J ∆ θ + ∇ x f ( x ∗ (ˆ θ ) , ˆ θ )∆ x (2.33) H ∆ x ≈ −J ∆ θ (2.34)where, (cid:101) x ∈ [ x ∗ (ˆ θ ) , x ∗ (ˆ θ ) + ∆ x ] (2.35)[ H ] ij = [ ∇ x f ( x ∗ (ˆ θ ) , ˆ θ )] ij (2.36)= ∂ [ ∇ x f ( x ∗ ; ˆ θ )] i ∂ x j (2.37)Now from [38], the gradient of the least squares function f can be written as ∇ x f ( x ∗ , ˆ θ ) = M T Φ( x ∗ , ˆ θ ) (2.38)[ M ] ij = ∂ Φ i ( x ∗ , ˆ θ ) ∂ x j (2.39)= ∂ F i ( x ∗ , ˆ θ ) ∂ x j if i (cid:54)∈ I ∂ψ ( x i , F i ( x ∗ , ˆ θ )) ∂ x j if i ∈ Z ∂ψ ( x i , F i ( x ∗ , ˆ θ )) ∂ x j otherwise (2.40)= ∂ F i ( x ∗ , ˆ θ ) ∂ x j if i (cid:54)∈ I i ∈ Z ∂ψ i ∂ x j otherwise (2.41)12hich is the form of M defined in algorithm 1. Also H = ∇ x f ( x ∗ ; ˆ θ ) = M T M + n (cid:88) i =1 Φ i ( x ∗ ; ˆ θ ) ∇ x Φ i ( x ∗ ; ˆ θ ) (2.42)= M T M (2.43)where the second term vanishes since we have from Theorem 5 that each termof Φ individually vanishes at the solution. Now J = ∇ x θ f ( x ∗ ; ˆ θ ) (2.44) J ij = ∂ [ ∇ x f ( x ∗ ; ˆ θ )] i ∂θ j (2.45)= ∂∂θ j (cid:32) n (cid:88) k =1 [ ∇ x Φ( x ∗ ; ˆ θ )] ki Φ k ( x ∗ ; ˆ θ ) (cid:33) (2.46)= n (cid:88) k =1 (cid:32) ∂ [ ∇ x Φ( x ∗ ; ˆ θ )] ki ∂θ j Φ k ( x ∗ ; ˆ θ ) + [ ∇ x Φ( x ∗ ; ˆ θ )] ki ∂ Φ k ( x ∗ ; ˆ θ ) ∂θ j (cid:33) (2.47)= n (cid:88) k =1 M ki N kj = M T N (2.48)where the first term vanished because Φ i are individually zeros, and we define N ij = ∂ Φ i ( x ∗ ; ˆ θ ) ∂θ j (2.49)= ∂ F i ∂ x j if i (cid:54)∈ I ψ ( x ∗ i ; F ∗ i ) ψ b ( x ∗ i ; F ∗ i ) ∂ F i ∂θ j if i ∈ Z ψ b ( x ∗ i ; F ∗ i ) ∂ F i ( x ∗ ;ˆ θ ) ∂θ j otherwise (2.50)= ∂ F i ∂θ j if i (cid:54)∈ I i ∈ Z ψ b ( x ∗ i ; F ∗ i ) ∂ F i ( x ∗ ;ˆ θ ) ∂θ j otherwise (2.51)which is the form of N defined in algorithm 1. By assumption 3, we have aunique minimum in the neighborhood of x ∗ where the gradient vanishes. So wehave from (2.34), (2.43) and (2.48) H ∆ x = −J ∆ θ (2.52) M T M ∆ x = −M T N ∆ θ (2.53)∆ x solves the above equation, if it solves M ∆ x = −N ∆ θ (2.54)13y defining T as the solution to the linear system of equations MT = N (2.55)∆ x = −T ∆ θ (2.56)and we have the above first-order approximation. From [41], we know that ifsome vector x has covariance C, then for a matrix A , the vector A x will havecovariance A C A T . So we have. C ov(∆ x ) ≈ T C ov(∆ θ ) T T (2.57)Thus we approximate the covariance of ∆ x in algorithm 1. (cid:4) The matrix M could potentially not have full rank, for example, when Z is non-empty, i.e., when we have weak complementarity terms. But we notethat M T M is the Hessian of f and the null-space of the Hessian correspondsto the directions where the gradient doesn’t change for small perturbations. Soin a first-order sence, small perturbations in those directions do not move thesolution, keeping the method robust even under weak complementarity. Thisis further confirmed by computational experiments detailed in Section 5.2 andAppendix D where we have cases with weak complementarity terms but nosignificant error.Further, when M is singular, we use Moore-Penrose pseudoinverse, M † tosolve the system of equations (2.55). Among possibly infinite solutions thatminimize the error (cid:107)MT − N (cid:107) , M † N gives the solution that minimizes (cid:107)T (cid:107) [5]. This would lead us to identifying the smallest step ∆ x that could be taken toreach the perturbed solution, up to first-order approximation and hence give themost conservative estimate of the uncertainty. Uniqueness and existence of M † is guaranteed and it can be computed efficiently. For computational purposesthe matrix T in the above equation has to be calculated only once, irrespectiveof the number of scenarios for which we would like to run for the covarianceof θ . Thus if x ∈ R n , θ ∈ R m and we want to test the output covariance for k different input covariance cases, the complexity is equal to that of solving asystem of n linear equations m times as in (2.55), and hence is O ( mn ). i.e., the complexity is quadratic in the number of output variables, linear in the14umber of input parameters and constant in the number of covariance scenarioswe would like to run.In Theorem 12 below, we prove that the error in the approximation of theo-rem 11 can be bounded using the condition number of the Hessian. We need thefollowing assumption that the condition number of the Hessian of f is boundedand the Hessian is Lipschitz continuous. Assumption 4.
At the known solution of the complementarity problem ofinterest ( θ = ˆ θ ),1. The condition number of the Hessian of f defined is finite and equals to κ H
2. The Hessian of f is Lipschitz continuous with a Lipschitz constant L ( x ∗ ; θ ). Theorem 12.
With assumption 4 holding, the error in the linear approximation2.34 for a perturbation of (cid:15) is o ( (cid:15) ) .Proof. Since ∇ f is Lipschitz continuous on both x and θ , we can write for (cid:101) x near x ∗ , (cid:13)(cid:13)(cid:13) ∇ x f ( x ∗ , ˆ θ ) − ∇ x f ( (cid:101) x , ˆ θ ) (cid:13)(cid:13)(cid:13) ≤ L ( x ∗ ; θ ) (cid:107) x ∗ − (cid:101) x (cid:107) (2.58) ≤ L ( x ∗ ; θ ) (cid:107) ∆ x (cid:107) (2.59) (cid:101) H = ∇ x f ( (cid:101) x , ˆ θ ) (2.60)= H + ε H (2.61)where (cid:107) ε H (cid:107) ≤ L ( x ∗ ; θ ) (cid:107) ∆ x (cid:107) . Applying the Lipschitz continuity on θ , (cid:13)(cid:13)(cid:13) ∇ θ ∇ x f ( x ∗ , (cid:101) θ ) − ∇ θ ∇ x f ( x ∗ , ˆ θ ) (cid:13)(cid:13)(cid:13) ≤ L ( x ∗ ; θ ) (cid:13)(cid:13)(cid:13)(cid:101) θ − ˆ θ (cid:13)(cid:13)(cid:13) (2.62) ≤ L ( x ∗ ; θ ) (cid:107) ∆ θ (cid:107) (2.63) (cid:101) J = ∇ θ ∇ x f ( x ∗ , (cid:101) θ ) (2.64)= J + ε J (2.65)where (cid:107) ε J (cid:107) ≤ L ( x ∗ ; θ ) (cid:107) ∆ θ (cid:107) . Thus the equation (cid:101) H ∆ x = (cid:101) J ∆ θ (2.66)is exact, even if we cannot compute (cid:101) H and (cid:101) J exactly. Now the error in inverting (cid:101) H is bounded by the condition number [26, Ch. 5]. (cid:13)(cid:13)(cid:13) H − − (cid:101) H − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H − (cid:13)(cid:13)(cid:13) ≤ κ H (cid:107) ε H (cid:107) (cid:107) (cid:101) H (cid:107) − κ H (cid:107) ε H (cid:107) (cid:107) (cid:101) H (cid:107) (2.67)15ssuming κ H (cid:107) ε H (cid:107) (cid:28) (cid:107)H(cid:107) , the above equation becomes (cid:13)(cid:13)(cid:13) H − − (cid:101) H − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H − (cid:13)(cid:13)(cid:13) ≤ κ H (cid:107) ε H (cid:107) (cid:13)(cid:13)(cid:13) (cid:101) H (cid:13)(cid:13)(cid:13) (2.68) ⇒ (cid:13)(cid:13)(cid:13) (cid:101) H − − H − (cid:13)(cid:13)(cid:13) ≤ κ H (cid:13)(cid:13)(cid:13) (cid:101) H − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H (cid:13)(cid:13)(cid:13) (cid:107) ε H (cid:107) (2.69) ⇒ (cid:13)(cid:13)(cid:13) (cid:101) H − (cid:101) J − H − J − H − ε J (cid:13)(cid:13)(cid:13) ≤ κ H (cid:13)(cid:13)(cid:13) (cid:101) H − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H (cid:13)(cid:13)(cid:13) ε H (cid:107)J (cid:107) + κ H (cid:13)(cid:13)(cid:13) (cid:101) H − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H (cid:13)(cid:13)(cid:13) (cid:107) ε H (cid:107) (cid:107) ε J (cid:107) (2.70) ⇒ (cid:13)(cid:13)(cid:13) (cid:101) H − (cid:101) J − H − J (cid:13)(cid:13)(cid:13) ≤ k (cid:107) ∆ x (cid:107) + k (cid:107) ∆ θ (cid:107) (2.71)with k = κ H L ( x ∗ ; θ ) (cid:13)(cid:13)(cid:13) (cid:101) H − (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:101) H (cid:13)(cid:13)(cid:13) (cid:107)J (cid:107) (2.72) k = L ( x ∗ ; θ ) (cid:13)(cid:13) H − (cid:13)(cid:13) (2.73)Thus we have from (2.71), that the error in the approximation done in algorithm1 is bounded. (cid:4)
3. Stochastic Sensitivity Analyses
In this section, we quantify the sensitivity of the variance of the solutionto variance in each of the input parameters. To achieve this, we define totallinear sensitivity and how it can be approximated using the matrix T derived in(2.55) . We then proceed to prove that these quantities also bound the maximumincrease in uncertainties of the output. Definition 13.
Given a function f : R m (cid:55)→ R n , the total linear sensitivity , β d ∈ R + of a dimension d ≤ m ; d ∈ N at a point x ∈ R m is defined for δ > β d = inf (cid:26) α : (cid:12)(cid:12)(cid:12)(cid:12) (cid:107) f ( x + δe d ) (cid:107) − (cid:107) f ( x ) (cid:107) (cid:12)(cid:12)(cid:12)(cid:12) ≤ δα + o (cid:0) δ (cid:1)(cid:27) (3.1)where e d is the d -th standard basis vector.This is a bound on the distance by which the function value can move for asmall perturbation in the input. Now we look at the solution to the parametrized16omplementarity problem in (2.1) as a function from the space of parametertuples to the space of solution tuples and bound the change in solution for asmall perturbation. The next proposition shows how the total linear sensitivitycan be calculated from the linear approximation matrix T derived earlier. Proposition 14.
Suppose we know, G ∈ R n × m such that G ij = ∂ f i ( x ) ∂ x j , then β d = (cid:112) ( (cid:80) ni =1 G id ) Proof.
See Appendix A (cid:4)
The above proposition proves that the T matrix obtained in (2.55) is suffi-cient to approximate the total linear sensitivity. The following result suggestshow the total linear sensitivity can approximate the total variance in the outputvariables. Theorem 15.
Given a function f : R m (cid:55)→ R n and β d , the increase in the totaluncertainty in the output, i.e., the sum of variances of the output variables, fora small increase of the variance of an input parameter, σ d of x d is approximatedby β d σ d .Proof. Let E d be the matrix of size m × m with zeros everywhere except the d -thdiagonal element, where it is 1. Given C = C ov( x ( ω )), for a small perturbation σ in the variance of x d , the covariance of f ( x ) changes as follows. C ∗ ≈ ∇ x f C ∇ x f T (3.2) C ∗ + ∆ C ∗ ≈ ∇ x f ( C + σ E d ) ∇ x f T (3.3)= C ∗ + σ ∇ x f E d ∇ x f T (3.4)[∆ C ∗ ] ij ≈ σ [ ∇ x f ] id [ ∇ x f ] jd (3.5) n (cid:88) i =1 [∆ C ∗ ] ii ≈ σ β d (3.6)which is the total increase in variance. The off-diagonal terms do not affectthe total uncertainty in the system because, the symmetric matrix C can bediagonalized as QDQ T , where Q is a rotation matrix, and the trace is invariantunder orthogonal transformations. (cid:4) With the above result, we can determine the contribution of each inputparameter to the total uncertainty in the output.17 . Application to optimization
To illustrate the application of this method explained in algorithm 1, weuse it to derive an approximation for the covariance of the solution of certaincanonical optimization problems. The goal of this section is to walk the readerthrough a simple application of the method to develop intuition of the analysisand results.To start with, we assume conditions on the differentiability and convexity ofthe objective function.
Assumption 5.
The objective function f ( x ; θ ) is strictly convex in x and istwice continuously differentiable in x and θ .In the theorem below, we approximate the covariance of the decision vari-ables of a convex optimization with uncertainties in the linear term and withonly linear equality constraints. Theorem 16.
With assumption 5 holding, the covariance of the primal anddual variables at the optimum of the problem,
Minimize x f ( x ; θ ) = g ( x ) + c ( θ ) T x (4.1) subject to Ax = b ( θ ) ( y ) (4.2) where θ = θ ( ω ) are random parameters with covariance C , is first-order approx-imated by T C T T where T = (cid:18) ∇ x g ( x ∗ ) A T A (cid:19) − (cid:18) −∇ θ c ( θ ) ∇ θ b ( θ ) (cid:19) (4.3) Proof.
For the given optimization problem, because of assumption 5 and linearindependence constraint qualification (LICQ), the KKT conditions are necessaryand sufficient for optimality. The KKT condition satisfied at a solution ( x ∗ , y ∗ )for the problem are given by ∇ x g ( x ∗ ) + c ( θ ) + A T y ∗ = 0 (4.4) A x ∗ = b ( θ ) (4.5)for some vector y so that the equation is well defined. Suppose from there, θ isperturbed by ∆ θ , we have ∇ x g ( x ∗ ) + c ( θ + ∆ θ ) + A T y ∗ ≈ ∇ θ c ( θ )∆ θ (4.6) A x ∗ − b ( θ + ∆ θ ) ≈ −∇ θ b ( θ )∆ θ (4.7)18ow we need to find ∆ x and ∆ y such that ∇ x g ( x ∗ + ∆ x ) + c ( θ + ∆ θ ) + A T ( y ∗ + ∆ y ) ≈ A ( x ∗ + ∆ x ) − b ( θ + ∆ θ ) ≈ ∇ x g ( x ∗ )∆ x + A T ∆ y ≈ ∇ θ c ( θ )∆ θ (4.10) A ∆ x ≈ −∇ θ b ( θ )∆ θ (4.11)The above conditions can be compactly represented as (cid:18) ∇ x g ( x ∗ ) A T A (cid:19) (cid:18) ∆ x ∆ y (cid:19) = (cid:18) ∇ θ c ( θ ) −∇ θ b ( θ ) (cid:19) ∆ θ (4.12)If A has full rank, then the above matrix is non-singular. So the change in thedecision variables x and the duals y can be written as a linear transformationof the perturbation in the random parameters. And we now have C ov (cid:18) ∆ x ∆ y (cid:19) = T C ov( θ ) T T (4.13) T = (cid:18) ∇ x g ( x ∗ ) A T A (cid:19) − (cid:18) −∇ θ c ( θ ) ∇ θ b ( θ ) (cid:19) (4.14) (cid:4) In the corollary below, we show that the method suggested is accurate ( i.e., has zero error) for an unconstrained quadratic optimization problem with un-certainty in the linear term.
Corollary 17.
For an optimization problem with uncertainty of objectives ofthe form, f ( x ; θ ) = 12 x T G x + θ ( ω ) T x (4.15) where G is positive definite, the approximation method has zero error. In otherwords, the obtained covariance matrix is exact.Proof. See Appendix A (cid:4)
5. Application to a general oligopoly market
We now present an example of a complementarity problem in a natural gasoligopoly and show how the methods developed in this paper can be applied.19 .1. Problem Formulation and results
Consider k producers competitively producing natural gas in a Nash-Cournotgame. Let the random unit costs of production be γ i ( ω ) , i ∈ { , . . . , k } . Also,let us assume that the consumer behavior is modeled by a linear demand curve P ( ˜ Q ) as follows. P = a ( ω ) + b ( ω ) ˜ Q (5.1)where P is the price the consumer is willing to pay, ˜ Q is the total quantity ofthe natural gas produced and random variables a ( ω ) > , b ( ω ) < ∀ ω ∈ Ω.Suppose the producers are maximizing their profits, then the Nash equilibriumcan be obtained by solving the following complementarity problem [12, 21].0 ≤ Q i ⊥ F i ( Q ) = γ i − a − b k (cid:88) j =1 Q k − bQ i ≥ a, b, γ i correspond to θ and Q i correspond to x in (2.1) with I = { , , . . . , k } . In the current numerical example, let us consider a duopolywhere k = 2. Let E (cid:16) γ γ a b (cid:17) T = (cid:16) − (cid:17) T (5.3)Solving the complementarity problem deterministically with the above parame-ter values, we get Q and Q to be 4 and 5 respectively. We use the C-function ψ min ( x, y ) = min( x, y ) for this example to get M = N = − −
130 1 − − (5.4)Now we have from (2.55) T = M − N = 13 − − − − − − (5.5)Having obtained T , we attempt to get insight on how uncertainties in variousinput parameters propagate through the model causing uncertainty in the equi-librium quantities. If we assume that all these parameters, viz. γ , γ , a, b have a200% coefficient of variation and are all uncorrelated, then the covariance matrixof the input is C = .
04 0 0 00 0 .
01 0 00 0 2 .
25 00 0 0 0 . (5.6)Then the covariance matrix of the solution would be C ∗ = T C T T = . . . . (5.7)The standard deviation of the produced quantities are 0.65(= √ . √ . C has the third and fourth diagonal term of C as zero. In such ascenario, we would expect the decrease in the quantity of production of oneplayer to cause an increase in the quantity of production of the other and viceversa, caused by re-adjustment of market share. We can see this effect by com-puting the covariance of the solution as T C T . The solution thus obtainedshows that the produced quantities are negatively correlated with a correla-tion of − −
62% correlation. The standard deviations in the pro-duced quantities have also dropped to about 2.9% and 1.2%. Thus we not only21btain insight about the uncertainties in the output, but also the correlation be-tween the output parameters. From an energy market policy maker’s perspectivethis is crucial information as it helps identifying the regions where increase ordecrease in production, consumption, price, pipeline flows and infrastructuralexpansions occur simultaneously and where they change asynchronously. Nowwe calculate the sensitivity of each of the input parameters to identify the pa-rameter that causes maximum uncertainty in the output. The values for β foreach of the four parameters γ , γ , a, b are calculated below. β = 13 (cid:16) √ √ √ √ (cid:17) T = (cid:16) .
745 0 .
745 0 .
471 6 . (cid:17) T (5.8)Thus we see that the solution is more sensitive to the slope of the demandcurve than to say production cost. Strictly speaking, if we define the variancein equilibrium as the sum of the variance of all output variables, this says, aunit increase in variance of the slope of the demand curve will be magnifiedabout 41 times (6 . ) variance in the equilibrium. However, a unit increase inthe variance of the production cost only increases the variance in equilibriumby 0.556 (0 . ) units. We used a Monte-Carlo based method as a comparison against our approxi-mation method to compute covariance in the decision variables. To achieve this,we modeled the oligopoly complementarity problem mentioned in (5.2) varyingthe number of players, and hence the number of random parameters and thedecision variables. For the Monte-Carlo simulation based approach, a symmet-rically balanced stratified design [43] is used with each dimension divided intotwo strata. With increasing number of random parameters and equilibriumvariables, Monte-Carlo methods become increasingly inefficient as the numberof simulations required grows exponentially. A comparison of the time taken inan processor to solve theabove oligopoly problem with varying number of players is shown in Fig. 2a.22 - - -
10 20 30 40 50
Number of producers C o m pu t a t i on t i m e i n m i nu t e s Approximation timeMonte Carlo timeMonte Carlo time(Extrapolated)
Computation time (a) Run time comparison (b) Error comparison
Figure 2: Computational experiments comparing Monte-Carlo methods and First-order ap-proximation methodFigure 3: Regional disaggregation of United States and Mexico. Source: [2] and U.S. EnergyInformation Administration
Despite developments in algorithms to solve complementarity problems, the saidexponential growth in the number of sample points required in a Monte-Carlobased approach deters the computational speed. A problem with as few as 25uncertain variables takes about 2 hours to solve and one with 30 uncertain vari-ables takes about seven days to solve using Monte-Carlo based approaches whileit takes few seconds to minutes in the first-order approximation method. Fig 2bcompares the error between 5 rounds of Monte-Carlo simulation and the first-order approximation method. More details on these computational experimentsare provided in Appendix D. 23 . Application to North American Natural Gas Market
In Mexico, motivation to move from coal to cleaner energy sources creates anincreasing trend in natural gas consumption, particularly in the power sector.Technology change, including fracking, has made natural gas available at alow cost. This resulted in increased production and higher proven reserves inthe U.S. Therefore, the US is expected to become a net exporter of NaturalGas (increasing pipelines exports to Mexico and LNG) during the next years[16, 18]. The North American Natural Gas Model (NANGAM) developed in[18] analyzes the impacts of cross border trade with Mexico. NANGAM modelsthe equilibrium under various scenarios by competitively maximizing the profitsof suppliers and pipeline operators, and the utility of consumers, resulting in acomplementarity problem. The model also uses the Golombek function [24, 27]to model the increase in marginal cost of production when producing close tocapacity. The formal description of the model is provided in Appendix B.For the model in this paper, which is motivated by NANGAM, we havedisaggregated the United States into 9 census regions (US1-9) and Alaska [2].Mexico is divided into 5 regions (MEX1-5). A map showing this regional disag-gregation is shown in Fig. 6. Further Canada is divided into two zones, CanadaEast (CAE) and Canada West (CAW). The model has 13 suppliers, 17 con-sumers, 17 nodes, and 7 time-steps. This amounts to 12,047 variables (primaland dual) and 2023 parameters. The gradient matrix of the complementarityfunction would contain 12 , elements and a Hessian matrix will have 12047 elements which is more than 1700 trillion floating point variables. We needefficient methods to handle these large objects. We observe, however, that thedependence of each component of the complementarity function is limited to fewvariables, thus making the gradient matrix sparse. Efficient sparse matrix toolsin scipy [31] are used along with a python class we specially built to handle asparse multi-dimensional array. The details of this class are given in AppendixC. This model is calibrated to match the region-wise production and consump-24ion data by adjusting the parameters of the demand curve, supply curve andthe transportation cost. The source for the projected numbers are the sameas the ones in Table 2 of [18]. The parameters of the demand curve were cho-sen in such a way that an elasticity of 0.29 is maintained at the solution to beconsistent with [15]. We used the method developed in algorithm 1 to understand the propaga-tion of uncertainty in the model. The covariance for each parameter acrossyears is obtained by fitting a Wiener process to the parameter value. This ischosen to mimic the Markovian and independent increment properties of marketparameters. Thus we have for any parameter dθ ( t ) = dµ θ ( t ) + σ θ dB ( t ) (6.1)where µ θ is calibrated, σ θ is chosen to be 1% of the average value of µ θ inthe analyzed period and B ( t ) is the standard Brownian motion. The diffusionparameter σ θ is assumed to be independent of time. Additionally to understandthe effect of greater uncertainty in US7, that accounts for about 40% of the totalproduction in the continent, the parameters of production cost are assumed tohave 5 times the variance than in any other region. The deterministic version of the problem is solved using the PATH algorithm[13] by assuming a mean value for all random parameters. Following this, algo-rithm 1 was applied and the T matrix defined in (2.55) is obtained by solvingthe linear system of equations using a Moore-Penrose pseudoinverse [26]. In thefollowing paragraph, we discuss some of the results obtained in this study.The heat map on Fig. 4a shows the coefficient of variation (standard devia-tion divided by mean) in consumer price in each year caused by the uncertaintyin parameters as mentioned in subsection 6.1. We notice that this large uncer-tainty in production costs of US7 caused relatively small uncertainties in theconsumer price. This is partially due to the availability of resources in US8 and25 a) Coefficient of variation in Price (b) Covariance of Produced quantity Figure 4: Covariance results (a) Price sensitivity to demand (b) Parameter sensitivity comparison
Figure 5: Sensitivity results
CAW to compensate for the large uncertainty in US7. The fact that it is actu-ally US8 and CAW that compensate for this uncertainty is known by looking atthe covariance plot on Fig. 4b which shows large correlation between US7 andUS8 and also between US7 and CAW.Fig. 5 shows the sensitivity of the solution to various input parameters.The graph on Fig. 5a shows the sum total change in uncertainty in price for a1% fluctuation in the demand curve of consumers. We notice that the price isparticularly sensitive to changes in demand in Mexico. This reflects the increas-ing concern about growing exports (both LNG and pipeline) that are likely toresult in higher consumer prices in the U.S. We also note that fluctuations indemand at nodes where production facilities are not available (MEX1, MEX3,MEX4) cause greater uncertainty in price. This is because, for regions with26 production facility in the same node, the production facility produces moreto cater the demand at that node and there is little effect in the flows and inthe prices at other nodes. This is also contingent to sufficient pipeline capac-ity. Larger changes in demand for regions with limited pupeline capacity (e.g.MEX1) may result in major changes in price. However a perturbation to thedemand at a node with no production unit causes the flows to alter to haveits demand catered. This affects natural gas availability elsewhere and causeslarger fluctuations in price. The tornado plot on Fig. 5b sorts the parametersin decreasing order of their effect on the uncertainty of the solution.The plot Fig. 5b shows the total change of the equilibrium if a parameter(e.g., the demand intercept) is shifted by 1% from its original value. Note thatthe results are plotted in logarithmic scale and are sorted in decreasing order. Ingeneral, our results suggest that parameters of consumers and producers playa major role on the equilibrium and hence a small perturbation have a largeeffect on the solution. In particular, the intercept and slope of the demandcurve and the linear cost parameter are the three most significant parameters.Interestingly, the expansion cost parameters (pipeline as well as production ex-pansion) have a lower effect on the solution equilibrium. As it was describedon Fig 5a, uncertainties in demand significantly affect regions with no or lowproduction capacities. Fig 5b corroborates that changes to the demand affectsthe equilibrium the most. The results also indicate that the infrastructure ex-pansion happens independently of changes in expansion cost. Natural gas pricespaid by consumers account for cost expansions and transportation. Therefore,the level of expansion is then driven by changes on demand. Hence, if policychanges need to be implemented in order to, for instance, to increase economicactivity or reduce carbon emissions, respectively subsidizes or taxes the down-stream or upstream ends of market rather than the mid-stream players to havelarger impacts. However, a policy maker who is interested in generating rev-enue without much impacts on the equilibrium should tax fuel transportationor infrastructure expansion for the greatest benefit.27 . Conclusion and Future work
In this paper, we developed a method to approximate the covariance of theoutput of a large-scale nonlinear complementarity problem with random inputparameters using first-order approximation methods. We extended this methodto general optimization problems with equality constraints. We then developedsensitivity metrics for each of the input parameters quantifying their contribu-tion to the uncertainty in the output. We used these tools to understand thecovariance in the equilibrium of the North American natural gas Market. Themethod gave insights into how production, consumption, pipeline flows, priceswould vary due to large uncertainties. While the variances identified the regionsthat are affected the most, the covariance gave information about whether thequantity will increase or decrease due to perturbation in the input. We alsoobtained results on the sensitivity of price uncertainty to demand uncertaintyin various nodes. We then quantified the contribution of each input parameterto the uncertainty in the output. This in turn, helps in identifying the regionsthat can have large impacts on equilibrium.We note that the method is particularly useful for large-scale nonlinear com-plementarity problems with a large number of uncertain parameters, which makeMonte-Carlo simulations intractable. It is robust in approximating the solutioncovariance for small uncertainty in the inputs. It is also good in quantifying thesensitivity of the output (and its variance) to the variance of input parameters.However since all the above are obtained as an approximation based on first-order metrics, there is a compromise in the accuracy if the variances of the inputare large. The method works the best for problems involving a large number ofdecision variables and random parameters with small variance.We foresee expanding this work by using progressively higher order termsof the Taylor series to capture the nonlinearities more efficiently. To ensurecomputational feasibility, this would typically require us to have stronger as-sumptions on the sparsity of the Hessian and the higher-order derivatives. Thiswill also require analysis and stronger assumptions about higher-order moments28f the random parameters.
8. Acknowledgements
The model in this article is based in part on the multi-fuel energy equilibriummodel MultiMod [28]. The MultiMod was developed by Dr. Daniel Huppmannat DIW Berlin as part of the RESOURCES project, in collaboration with Dr.Ruud Egging (NTNU, Trondheim), Dr. Franziska Holz (DIW Berlin) and others(see http://diw.de/multimod ). We are grateful to the original developersof MultiMod for sharing the mathematical implementation, which we furtherextended as part of this work.The authors would also like to thank Dr. Donniell Fishkind, Departmentof Applied Mathematics and Statistics, Johns Hopkins University, Dr. MichaelFerris, Department of Computer Sciences, University of Wisconsin-Madison andthe participants of TAI conference 2016, MOPTA 2016 and INFORMS 2016 fortheir valuable comments and discussions.The authors would also like to thank the two anonymous reviewers whosecomments and suggestions improved this paper.
9. Proofs to certain lemmas and propositions
Proof of Lemma 3.
To show this, we first prove that every element in K (cid:48) indeedis in K ∗ . And then we prove for every element x (cid:54)∈ K (cid:48) , there exists some v ∈ K such that v T x < x in K (cid:48) . v T x = n (cid:88) i =1 v i x i = (cid:88) i ∈I v i x i + (cid:88) i (cid:54)∈I v i x i ≥ K (cid:48) ⊆ K ∗ .Now to show the reverse containment, suppose there is x ∈ R n ; x (cid:54)∈ K (cid:48) . Thismeans, we either have1. at least one index j ∈ I such that x j < j (cid:54)∈ I such that x j (cid:54) = 029ow, v T x = n (cid:88) i =1 v i x i = v j x j + (cid:88) i (cid:54) = j v i x i (9.2)In the first case, choose v such that [ v ] i = 0 for i (cid:54) = j and v j = 1. Clearly v ∈ K and for this choice of v , the above sum is negative, showing x (cid:54)∈ K ∗ . In thesecond case, choose v such that [ v ] i = 0 for i (cid:54) = j and v j = − sgn( x j ). Clearly v ∈ K and for this choice of v , the above sum is negative, showing x (cid:54)∈ K ∗ . Thuswe show ( K (cid:48) ) c ⊆ ( K ∗ ) c , which implies the reverse containment and completesthe proof. (cid:4) Proof of Proposition 5.
Since x ∗ ≡ x ∗ ( θ ) solves the problem, following from therequirement that F ( x ∗ ) ∈ K ∗ and lemma 3, if i (cid:54)∈ I , F i ( x ∗ ; θ ) = 0.For i ∈ I , x ∗ ∈ K ⇒ x ∗ i ≥ F ( x ∗ ) ∈ K ∗ ⇒ F i ( x ∗ ) ≥
0. Also from therequirement x ∗ T F ( x ∗ ) = 0, one of the above two quantities should vanish foreach i ∈ I . But C-functions are precisely functions that vanish when both theirarguments are non-negative and of them equal zero. So ψ i ( x ∗ , F i ( x ∗ )) = 0.Thus each coordinate of Φ is individually zero, which makes f ( x ∗ ) vanish, whichis the smallest value f can take. Thus x ∗ is a global minimum of f . (cid:4) Proof of Proposition 6.
Since a solution exists for the NCP, we know by propo-sition 5 that the minimum value f can take is 0. Suppose we have x ∗ ∈ R n suchthat f ( x ∗ ; θ ) = 0. Since f is sum of squares, this can happen only if each of theindividual terms are zero. This means for i (cid:54)∈ I , F i ( x ∗ ) = 0.Now since ψ i ( x i , F i ( x ∗ )) = 0 for i ∈ I , we know F i ( x ∗ ) ≥
0. This combinedwith the previous point implies F ( x ∗ ) ∈ K ∗ .Also from the fact that ψ i ( x ∗ i ; F i ( x ∗ ; θ )) = 0 for i ∈ I , we know that x ∗ ∈ K .It also implies that x ∗ i F i ( x ∗ ) = 0 for i ∈ I . Thus x ∗ T F ( x ∗ ; θ ) = n (cid:88) i =1 x ∗ i F i (9.3)= (cid:88) i ∈I x ∗ i F i + (cid:88) i (cid:54)∈I x ∗ i F i (9.4)= 0 + 0 = 0 (9.5)This implies x ∗ ( θ ) ⊥ F ( x ∗ ; θ ) and x ∗ ( θ ) solves the complementarity problem. (cid:4) Proof of Corollary 8.
The set O = { ( a, b ) : ψ min ( a, b ) = 0 } \ { (0 , } is the pos-itive coordinate axes except the origin. Rewriting this C-function as ψ min ( a, b ) = (cid:26) a if a < bb otherwise (9.6)30e see that the second derivative of ψ min vanishes at O . Also we observe thatlim ( a,b ) → (0 , ψ ( a, b ) ∂ ψ ( a, b ) ∂a∂b = 0 (9.7)since the second derivative is 0 everywhere except along the line a = b . (cid:4) Proof of Corollary 9.
For ψ = ψ F B , we have assumption 1 of proposition 7satisfied by [17]. For assumption 2,lim ( a,b ) → (0 , ψ ( a, b ) ∂ ψ ( a, b ) ∂a∂b = lim ( a,b ) → (0 , (cid:16)(cid:112) a + b − a − b (cid:17) ab (cid:0) √ a + b (cid:1) (9.8)= 0 (9.9)Thus f is twice continuously differentiable at its zeros. The twice continuousdifferentiability elsewhere follows directly from the fact that ψ F B is twice con-tinuously differentiable everywhere except at the origin. This ensures that allthe terms in the derivative of the sum of squares exist and are finite. (cid:4)
Proof of Proposition 14.
By definition, for some admissible d , f ( x + δe d ) = f ( x ) + δGe d + o ( δ ) (9.10) ⇒ [ f ( x + δe d )] i = [ f ( x )] i + δG id + o ( δ ) (9.11) (cid:107) f ( x + δe d ) (cid:107) ≤ (cid:107) f ( x ) (cid:107) + (cid:107) δG .d (cid:107) + (cid:107) o ( δ ) (cid:107) (9.12)= (cid:107) f ( x ) (cid:107) + δ (cid:118)(cid:117)(cid:117)(cid:116)(cid:32) n (cid:88) i =1 G id (cid:33) + o ( δ ) (9.13)where G .d is the d -th column of G . Also we have from (9.11) for sufficientlysmall δ , (cid:107) f ( x + δe d ) (cid:107) ≥ (cid:107) f ( x ) (cid:107) − (cid:107) δG .d (cid:107) + (cid:107) o ( δ ) (cid:107) (9.14)= (cid:107) f ( x ) (cid:107) − δ (cid:118)(cid:117)(cid:117)(cid:116)(cid:32) n (cid:88) i =1 G id (cid:33) + o ( δ ) (9.15) (cid:4) Proof of Corollary 17.
For the problem to be well-defined, let G ∈ R n × n and θ ∈ R n . This makes ∇ x θ f ( x ; θ ) ∈ R n × n . ∇ x f ( x ; θ ) = G x + θ ( ω ) (9.16) ∇ x f ( x ; θ ) = G (9.17)[ ∇ x θ f ( x ; θ )] ij = I (9.18)31 igure 6: Regional disaggregation of United States and Mexico. Source: [2] and U.S. EnergyInformation Administration Due to absence of terms dependent on x in the last two equations, we have anexact equation, G ∆ x = ∆ θ (9.19)Due to the exactness of the above equation, we have T = G − (9.20) C ov(∆ x ) = T C ov(∆ θ ) T T (9.21)with no error. (cid:4)
10. Natural Gas Market - Complementarity Formulation
In this formulation, we assume we have a set of suppliers P, consumers Cand a pipeline operator. The players are located in a set of nodes N, and someof them are connected by pipelines A.Let also say that P n ⊆ P, C n ⊆ C are located in node n ∈ N. Let A n be thepipelines connected to node n . The symbols used here are explained in Table 1,2 and 3. Most of the analysis closely follow [18] and [14]. Random parametersare denoted by an ( ω ) beside them. The implementation of this problem is madeavailable in https://github.com/ssriram1992/Stoch Aprx cov .32 able 1: Sets Set Explanation Set Explanation
P Set of suppliers A Set of pipeline connections(arcs)C Set of consumers A + n Set of arcs from node n on which natural gas flows outN Set of nodes A − n Set of arcs from node n on which natural gas flows inY Set of periods Maximize (cid:88) Y df y ( ω ) (cid:40)(cid:88) C Q C pcny π cy − Gol (cid:0) Q P pny , CAP P py (cid:1) − π XP py ( ω )X P py − (cid:88) A + n π A ay Q A pay (10.1)subject toQ C pcny , Q P pny , Q A pay ≥ P py , CAP P py ≥ P pny ≤ α P CAP P py (cid:0) δ py (cid:1) (10.2a)CAP P py = ˆQ p + y (cid:88) i =1 X P pi (cid:0) δ py (cid:1) (10.2b) (cid:88) C n Q C pcny + (cid:88) A + n Q A pay = Q P pny (1 − L P py ( ω ))+ (cid:88) A − n Q A pay (1 − L A ay ( ω )) (cid:0) δ pny (cid:1) (10.2c)where Gol( . ) = ( l P py ( ω ) + g P py ( ω ))Q P pny + q P py ( ω )Q P pny + g P py ( ω )(CAP P py − Q P pny ) log (cid:32) − Q P pny CAP P py (cid:33) (10.3) Maximize (cid:88) Y df y ( ω ) (cid:40)(cid:88) A Q A ay (cid:0) π A ay − γ A ya ( ω ) (cid:1) − π XA ay ( ω )X A ay (cid:41) (10.4)33 able 2: Symbols - Variables Symbol ExplanationQuantities Q C pcny Quantity produced by p in n to send to c in year y Q P pny Total quantity produced by p in year y Q A pay Total quantity p choses to send by arc a in year y Q A ay Total quantity sent by a during year y Prices π cy Unit price paid by consumer C in year Yπ A ay Unit price of sending natural gas through a during year y Capacity X P py Production expansion in year y for supplier p X A ay Transportation capacity expansion in year y for arc a CAP P py Production capacity for supplier p in year y CAP A ay Transportation capacity for arc a in year y Table 3: Symbols - Parameters
Symbol ExplanationQuantities ˆQ p Initial capacity of production for supplier p ˆQ a Initial capacity of transportation for pipeline a Prices π XP py ( ω ) Price of capacity expansion for supplier pπ XA ay ( ω ) Price of capacity expansion for transportation arc a Losses L P py ( ω ) Percentage loss in production by supplier p in year yL A ay ( ω ) Percentage loss in transportation via arc a in year yα P Availability fraction of the production capacity
Consumer E C cy ( ω ) Intercept of the demand curve for consumer c in year yD C cy ( ω ) Slope of the demand curve for consumer c in year y df y ( ω ) Discount Factor for year y A ay , X A ay , CAP A ay ≥ A ay ≤ CAP A ay (cid:0) δ ay (cid:1) (10.5a)CAP A ay = ˆQ a + y (cid:88) i =1 X A ai (cid:0) δ ay (cid:1) (10.5b) π cy = E C cy ( ω ) + D C cy ( ω ) (cid:88) P Q C pcny ( π cy ) (10.6)It can be shown that the above said optimization problems are all convex withnon-empty interior. Hence the Karush-Kuhn Tucker conditions (KKT condi-tions) are necessary and sufficient for optimality. The KKT conditions arepresented below and they form the equations for the complementarity problemalong with the constraints above. − df y ( ω ) π cy + δ pny ≥ (cid:0) Q C pcny (cid:1) (10.7a) df y ( ω ) π XP py ( ω ) − y (cid:88) i =1 δ pi ≥ (cid:0) X P py (cid:1) (10.7b) df y ( ω ) π A ay + (cid:16) I a ∈ A + n − I a ∈ A − n (1 − L A ay ( ω )) (cid:17) δ pny ≥ (cid:0) Q A pay (cid:1) (10.7c) df y ( ω ) ∂ Gol ∂ Q P pny + δ py − δ pny (1 − L P py ( ω )) ≥ (cid:0) Q P pny (cid:1) (10.7d) df y ( ω ) ∂ Gol ∂ CAP P py + α P δ py − δ py ≥ (cid:0) CAP P py (cid:1) (10.7e)35 − df y ( ω ) π A ay + γ A ya ( ω ) + δ ay ≥ (cid:0) Q A ay (cid:1) (10.8a) df y ( ω ) π XA ay ( ω ) − y (cid:88) i =1 δ ai ≥ (cid:0) X A ay (cid:1) (10.8b) δ ay − δ ay ≥ (cid:0) CAP A ay (cid:1) (10.8c) Q A ay = (cid:88) P Q A pay ( π A ay ) (10.9)
11. N-dimensional Sparse array implementation
A general purpose Python class has been implemented to handle a sparse ndarray object. The class is a generalization of the scipy class coo matrix whichstores the array coordinates of each non-zero element in the array. We nowdescribe the details of the implementation. A continuously updated version ofthe class can be found at https://github.com/ssriram1992/ndsparse . The n-dimensional sparse array ( coo array ) can be initialized by any of thefollowing methods. • A tuple , which initializes the sparse array of the shape mentioned in the tuple and with zeros everywhere. • A dense ndarray which will be converted and stored as a coo array . • A matrix of positions and a 1 dimensional array of values where thematrix contains the positions of the non-zero elements and the vectorcontaining the non-zero values of those positions. In this case the shapeof the coo array would be the smallest ndarray that can store all theelements given. Optionally a tuple containing the shape of the ndarray can be given explicitly. 36 Another coo array whose copy is to be created.
The following methods and attributes are available in the coo array . • print(coo array) will result in printing the location of each of the non-zero elements of the array and their values. • coo array.flush(tol = 1e-5) will result in freeing the space used instoring any zero-elements or elements lesser than the tolerance, tol . Suchnumbers typically arise out arithmetic operations on coo array or poorinitialization. • coo array.size() returns the number of non-zero elements in the coo array . • coo array.shape returns the shape of the underlying dense matrix. • coo array.add entry(posn,val) and coo array.set entry(posn,val) both add a new non-zero element with the given value at the given po-sition. The difference however is that set entry() checks if a non-zerovalue already exists at the mentioned position, and if yes, overwrites it.This search makes set entry() slower compared to add entry() whichassumes that the previous value or the position is zero. Thus add entry() could potentially cause duplicates and ambiguity, if an illegal input isgiven. However in case the input is ensured to be legal, add entry() ismuch faster. • coo array.get entry(posn) returns the value at the given position. • coo array.swapaxes(axis1,axis2) is a higher dimensional generaliza-tion of matrix transposes where the dimensions that have to swapped canbe chosen. • coo array.remove duplicate at(posn,func=0) checks if there are mul-tiple values defined for a single position in the sparse array. If yes, theyare replaced by a single entry containing the scalar valued defined by func
37r passes them to a function defined in func and stores the returned value.Passing a function for the argument func is incredibly useful in performingarithmetic operations on coo array . • coo array.todense() returns a dense version of the coo array . • coo array.iterate() returns an iterable over the non-zero positions andvalues in the coo array .The above class is used extensively to handle high-dimensional sparse arrays re-sulting out of variables containing pipelines, viz., Q A pay , Q A ay , ˆQ p , ˆQ a , π A ay , X A ay , δ ay , δ ay and parameters with pipelines, viz., CAP A ay , π XA ay ( ω ) , L A ay ( ω ) , γ A ya ( ω ).
12. Computational Experiments
The single-node single-product oligopoly mentioned in Section 5.1 is usedfor the computational experiments. Experiments were done with varying thenumber of players and to avoid any advantage due to symmetry, a differentvalue was used for the cost of production for each of the n players. The valuesused for the computational study are given below. a = 500 (12.1) b = − . E [ c i ] = 100 + 3 i ( i = 1 , , . . . , n ) (12.3) V ar ( c i ) = 1 ( i = 1 , , . . . , n ) (12.4)For the time tests, the comparison was between a Monte-Carlo simulation in-volving 0 . × n random samples and the approximation method.For the error comparison, the Monte-Carlo simulation was done using max(100 , . × n ) samples. The process was repeated five times to show the relative differ-ences between Monte-Carlo solutions purely due to the randomness in sampling.Having obtained the covariance matrix of the solution, both by Monte-Carlosimulation as well as the first-order approximation, we compute the trace of the38 igure 7: Error comparison covariance matrix. As said in Section 3, this trace corresponds to the invarianttotal uncertainty in the solution. Each cycle of the Monte-Carlo simulation givesone value of the said quantity. Each of these points is shown as a blue dot inthe figure. The mean of all these values is considered as the zero line. The reddot corresponds to the value obtained by the approximation method. We hencedemonstrate the relative accuracy of our method with respect to Monte-Carlosimulations. We also note that the test involves complementarity problems withexclusively strong complementarity terms which happens when the number ofplayers is at most 15, and problems with both strong and weak complementarityterms, which happens when the number of players exceed 15. In either case, theerror is comparable to that obtained via Monte-Carlo simulation.For the purposes of this experiment, the first-order approximation methodand the sampling process were implemented in Python 2.7 while the complemen-tarity problems were solved using PATH algorithm accessed through Python-GAMS api. 39 eferencesReferences [1] Abada, I., Gabriel, S., Briat, V., and Massol, O. (2013). A generalized nash–cournot model for the northwestern european natural gas markets with a fuelsubstitution demand function: The gammes model. Networks and SpatialEconomics , 13(1):1–42.[2] AEO, A. E. O. (2015). US Energy Information Administration, 2015.[3] Agrawal, S., Ding, Y., Saberi, A., and Ye, Y. (2012). Price of Correlationsin Stochastic Optimization.
Operations Research , 60(1):150–162.[4] Anderson, E. J. and Xu, H. (2004). Nash equilibria in electricity mar-kets with discrete prices.
Mathematical Methods of Operations Research ,60(2):215–238.[5] Ben-Israel, A. and Greville, T. N. (2003).
Generalized inverses: theory andapplications , volume 15. Springer Science & Business Media.[6] Benedetti-Cecchi, L. (2003). The importance of the variance around themean effect size of ecological processes.
Ecology , 84(9):2335–2346.[7] Castillo, E., Conejo, A., Castillo, C., M´ınguez, R., and Ortigosa, D. (2006).Perturbation approach to sensitivity analysis in mathematical programming.
Journal of Optimization Theory and Applications , 128(1):49–74.[8] Castillo, E., Conejo, A. J., and Aranda, E. (2008). Sensitivity analysis incalculus of variations. some applications.
SIAM review , 50(2):294–312.[9] Chen, X. and Fukushima, M. (2005). Expected Residual MinimizationMethod for Stochastic Linear Complementarity Problems.
Mathematics ofOperations Research , 30(4):1022–1038.[10] Chen, Y., Cowling, P., Polack, F., Remde, S., and Mourdjis, P. (2017).Dynamic optimisation of preventative and corrective maintenance schedules40or a large scale urban drainage system.
European Journal of OperationalResearch , 257(2):494–510.[11] Christensen, A. and Siddiqui, S. (2015). A mixed complementarity modelfor the us biofuel market with federal policy interventions.
Biofuels, Bioprod-ucts and Biorefining , 9(4):397–411.[12] Cottle, R. W., Pang, J. S., and Stone, R. E. (2009).
The linear comple-mentarity problem , volume 60. Siam.[13] Dirkse, S. P. and Ferris, M. C. (1995). The path solver: a nommono-tone stabilization scheme for mixed complementarity problems.
OptimizationMethods and Software , 5(2):123–156.[14] Egging, R., Pichler, A., Kalvø, Ø. I., and WalleHansen, T. M. (2016). RiskAversion in Imperfect Natural Gas Markets.
European Journal of OperationalResearch .[15] EIA, E. I. A. (2012). Fuel Competition in Power Generation and Elasticitiesof Substitution.[16] EIA, E. I. A. (2016). U.s. natural gas exports to mexico continue to grow.[17] Facchinei, F. and Pang, J. S. (2007).
Finite-dimensional variational in-equalities and complementarity problems . Springer Science & Business Media.[18] Feijoo, F., Huppmann, D., Sakiyama, L., and Siddiqui, S. (2016). NorthAmerican natural gas model: Impact of cross-border trade with Mexico.
En-ergy , 112:1084–1095.[19] Ferris, M. C. and Pang, J. S. (1997). Engineering and Economic Applica-tions of Complementarity Problems.
SIAM Review , 39(4):669–713.[20] Fiacco, A. V. (2009). Sensitivity and stability in NLP: Approximation Sen-sitivity and Stability in NLP: Approximation.
Encyclopedia of Optimization ,pages 3454–3467. 4121] Gabriel, S. A., Conejo, A. J., Fuller, J. D., Hobbs, B. F., and Ruiz, C.(2012).
Complementarity modeling in energy markets , volume 180. SpringerScience & Business Media.[22] Gabriel, S. A., Kydes, A. S., and Whitman, P. (2001). The National En-ergy Modeling System: A Large-Scale Energy-Economic Equilibrium Model.
Operations Research , 49(1):14–25.[23] Gabriel, S. A., Zhuang, J., and Egging, R. (2009). Solving stochastic com-plementarity problems in energy market modeling using scenario reduction.
European Journal of Operational Research , 197(3):1028–1040.[24] Golombek, R., Gjelsvik, E., and Rosendahl, K. E. (1995). Effects of lib-eralizing the natural gas markets in Western Europe.
The Energy Journal ,pages 85–111.[25] G¨urkan, G. and Robinson, S. M. (1999). Sample-path solution of stochasticvariational inequalities.
Mathematical Programming, Series B , 84(2):313–333.[26] Horn, R. A. and Johnson, C. R. (2012).
Matrix analysis . Cambridgeuniversity press.[27] Huppmann, D. (2013). Endogenous production capacity investment innatural gas market equilibrium models.
European Journal of OperationalResearch , 231(2):503–506.[28] Huppmann, D. and Egging, R. (2014). Market power, fuel substitutionand infrastructure–a large-scale equilibrium model of global energy markets.
Energy , 75:483–500.[29] Hyett, K. L., Podosky, M., Santamaria, N., and Ham, J. C. (2007). Valuingvariance: the importance of variance analysis in clinical pathways utilisation.
Australian Health Review , 31(4):565–570.[30] Jiang, H. and Xu, H. (2008). Stochastic approximation approaches to thestochastic variational inequality problem.
IEEE Transactions on AutomaticControl , 53(6):1462–1475. 4231] Jones, E., Oliphant, T., Peterson, P., and Others (2015). { SciPy } : Opensource scientific tools for { Python } .[32] Kopanos, G. M., M´endez, C. A., and Puigjaner, L. (2010). Mip-baseddecomposition strategies for large-scale scheduling problems in multiproductmultistage batch plants: A benchmark scheduling problem of the pharmaceu-tical industry. European journal of operational research , 207(2):644–655.[33] Krantz, S. G. and Parks, H. R. (2012).
The implicit function theorem:history, theory, and applications . Springer Science & Business Media.[34] Kwak, B. M. and Lee, S. S. (1988). A complementarity problem formulationfor two-dimensional frictional contact problems.
Computers & structures ,28(4):469–480.[35] Lamm, M., Lu, S., and Budhiraja, A. (2016). Individual confidence in-tervals for solutions to expected value formulations of stochastic variationalinequalities.
Mathematical Programming , pages 1–46.[36] Luo, J., Hong, L. J., Nelson, B. L., and Wu, Y. (2015). Fully Sequen-tial Procedures for Large-Scale Ranking-and-Selection Problems in ParallelComputing Environments.
Operations Research , 63(5):1177–1194.[37] Mart´ın, S., Smeers, Y., and Aguado, J. A. (2015). A Stochastic TwoSettlement Equilibrium Model for Electricity Markets With Wind Generation.
IEEE Transactions on Power Systems , 30(1):1–13.[38] Nocedal, J. and Wright, S. (2006).
Numerical optimization . Springer Sci-ence & Business Media.[39] Ohno, K., Boh, T., Nakade, K., and Tamura, T. (2016). New approximatedynamic programming algorithms for large-scale undiscounted markov deci-sion processes and their application to optimize a production and distributionsystem.
European Journal of Operational Research , 249(1):22–31.4340] Oke, O., Huppmann, D., Marshall, M., Poulton, R., and Siddiqui, S. (2016).Mitigating environmental and public-safety risks of united states crude-by-railtransport.[41] Seber, G. A. F. and Lee, A. J. (2012).
Linear regression analysis , volume936. John Wiley & Sons.[42] Shanbhag, U. V. (2013).
Stochastic Variational Inequality Problems : Ap-plications , Analysis , and Algorithms . INFORMS.[43] Shields, M. D., Teferra, K., Hapij, A., and Daddazio, R. P. (2015). Refinedstratified sampling for efficient monte carlo based uncertainty quantification.
Reliability Engineering & System Safety , 142:310–325.[44] Siddiqui, S. and Christensen, A. (2016). Determining energy and climatemarket policy using multiobjective programs with equilibrium constraints.
Energy , 94:316–325.[45] Yumashev, D. and Johnson, P. (2017). Flexible decision making in the wakeof large scale nuclear emergencies: Long-term response.