© 2008 Ankur A. Kulkarni.
TOPICS IN STOCHASTIC OPTIMIZATION AND EQUILIBRIUM PROBLEMS
BY ANKUR A. KULKARNI B.Tech. Indian Institute of Technology, Bombay, 2006.
THESIS Submitted in partial fulfillment of the requirements for the degree of Master of Science in General Engineering in the Graduate College of the University of Illinois at Urbana-Champaign, 2008
Urbana, Illinois
Adviser: Asst. Prof. Vinayak “Uday” V. Shanbhag
Abstract
Stochastic programming forms a class of optimization problems that has seen much attention in the last decade. As has the class of equilibrium problems which have resulted largely from the analytic treatment of Nash games. This thesis will delve into some challenging areas of stochastic optimization and stochastic equilibrium problems. Our interest lies chiefly in algorithms for solving such problems. This thesis comprises of the following studies. Decomposition methods for stochastic nonlinear programs: In this paper, we consider a two-period stochastic nonlinear program and consider the use of sequential quadratic programming for the solution of such problems. Expectedly, the quadratic programming approximation may be arbitrarily large but has an attractive block-angular structure. To address such problems, we may employ between one of two QP solvers at our disposal. The first is based on the inexact cut version of the L-shaped method for stochastic quadratic programming [81, 82]. The second is a QP extension of the trust-region method suggested by Linderoth and Wright [56]. We also use a line-search in line with conventional SQP ideas. We compare the performance of the algorithms on a quadratic programming test problem set. Splitting method for Nash-Cournot equilibria under uncertainty: Inspired by a Cournot bidding model for electricity markets subject to uncertainty, we first show existence of Nash equilibria to a wide class of stochastic quadratic games. We also develop a splitting method that can be decomposed scenario-wise to compute such equilibria. We also indicate how our splitting method can serve as a subproblem solver for nonlinear stochastic Nash-Cournot games.
ii
To Aai, Pappa, Apoorva, and IIT Bombay.
iii
Acknowledgments
When I joined the University of Illinois I was unaware of what applied mathematical research entails. As a naive undergraduate in IIT Bombay I was fortunate to have Uday take me under his wing and be my adviser. He has been a guru in the true sense of the Sanskrit word: knowledgeable, open minded, motivating, and caring. He has great depth of knowledge in his area and is refreshingly open about areas of research. Whenever I have needed his help, even on the tiniest of things, he has always been accessible, patient and eager to help. It has been only two years, but I can already recall numerous energetic, inspiring discussions we have had. His contribution to my wellbeing has extended outside the academic and professional environment. A year ago when I was down with flu and entire campus was under a foot of snow, Uday was kind enough to drive me to and fro to the hospital. Every time we have worked late at night, he has been kind enough to offer me a ride home. It is these personal interactions and the care he has shown for my growth that has made working with him special. His understanding of my problems and patience in dealing with them has seen my naivety reduce significantly and has readied me for the Ph.D. ahead. The University of Illinois brought me a windfall by the name Anupama, my soulmate, and to be fiancee. She has been extremely supportive and incredibly patient, forgiving and understanding through the many hours of togetherness that this thesis has stolen from her. Even without her presence, this thesis would have been written, the results would have been published, and my degree would have been obtained. But the spirit that this work carries with itself would not have existed. This thesis has life and my life has beauty and optimism only because of her. Special thanks are due to roomie Pritam: those sambars and rassams shall remain an enduring memory; the IESE office staff: in particular Donna for all the effort she has put in for this thesis and for uncountable little things that I need her help with. Lastly, I thank my family, and IIT Bombay, to whom this thesis is dedicated.
iv
No words are enough to describe the love my family has showered on me. None of the above acknowledgments would have occurred without them. My teachers and friends at IIT Bombay taught me so many lessons which have shaped the way I think and the way I approach life. Thank you for everything.
v
Table of Contents
1 A Benders-SQP Framework for Nonlinear Programming under Uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Financial Planning . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Newsvendor Problem . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Stochastic Equilibrium Problems . . . . . . . . . . . . . . . . 1.2.4 Summary of Literature . . . . . . . . . . . . . . . . . . . . . . 1.2.5 Random Recourse . . . . . . . . . . . . . . . . . . . . . . . . 1.3 An SQP Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Sequential Quadratic Programming . . . . . . . . . . . . . . 1.3.2 A Linesearch SQP Method . . . . . . . . . . . . . . . . . . . . 1.3.3 Constructing QP Approximations . . . . . . . . . . . . . . . 1.3.4 A Parallelizable Linesearch . . . . . . . . . . . . . . . . . . . 1.3.5 A Sparse Quasi-Newton Update . . . . . . . . . . . . . . . . 1.4 Dual Decomposition Methods for the QP Approximations . . . . 1.4.1 L-shaped Method with Inexact Cuts . . . . . . . . . . . . . . 1.4.2 Trust-region L-shaped Method for Stochastic QPs . . . . . . 1.4.3 Recovering the Multipliers . . . . . . . . . . . . . . . . . . . . 1.5 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Scalability of the Algorithm . . . . . . . . . . . . . . . . . . . 1.5.2 Modifying the TR Update and Inexactness Update . . . . . 1.5.3 Comparing ILS and TR . . . . . . . . . . . . . . . . . . . . . . 1.6 Summary, Concluding Remarks and Future Directions . . . . . . . 1.7 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1 1 4 4 4 5 6 8 15 16 17 20 22 24 29 31 33 45 48 48 50 53 54 55
2 Stochastic Nash-Cournot Equilibrium Problems: Theoretical Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.2 Motivation and Background . . . . . . . . . . . . . . . . . . . . . . . 59 2.2.1 Source Problems . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.2.2 Outline of Literature . . . . . . . . . . . . . . . . . . . . . . . 61 2.3 The Deterministic Nash-Cournot Problem . . . . . . . . . . . . . . 71 2.4 The Stochastic Nash-Cournot Equilibrium Problem . . . . . . . . . 76 2.4.1 Formulating the Stochastic Nash-Cournot Problem . . . . . 77
vi
3 An Inexact-Newton Splitting (INS) Method . . . . . . . . . . 3.1 An Inexact Newton Framework . . . . . . . . . . . . . . . 3.2 A Matrix-Splitting Subproblem Solver . . . . . . . . . . . 3.2.1 Linear constraints and contraction mappings . . 3.3 Convergence Properties of the Algorithm . . . . . . . . . 3.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Scalability . . . . . . . . . . . . . . . . . . . . . . . 3.5 Summary and Conclusions . . . . . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. 82 . 83 . 87 . 90 . 93 . 93 . 94 . 96
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
98
Author’s Biography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
vii
1
A Benders-SQP Framework for Nonlinear Programming under Uncertainty ———————————–
1.1 Introduction In this paper1 , we focus on stochastic nonlinear programs. As the name suggests, these are optimization problems in which there is an element of stochasticity or uncertainty manifesting itself in the form of random parameters. We consider the problem (SNLP): SNLP
£ ¤ minimize f (x) + IEω h ω (y ω ) x,y
u(x) = 0 subject to
a ω (x) + d ω (y ω ) = 0 x, y ω ≥ 0,
∀ω∈Ω
where f , h ω , u, a ω and d ω are assumed to be continuously differentiable nonlinear functions for all ω ∈ Ω. The expectation IE(·) is taken with respect to a measure IP for the discrete random variable ω which belongs to Ω, a sample space of finite cardinality. Stochastic programming formulations were separately suggested by Dantzig [21] and Beale [7] for planning under uncertainty. The challenge in such problems lies in coping with the size of the problem, given that the sample space Ω may be become arbitrarily large. Hence computational methods have centered around decomposing the problem into a set of tractable (smaller) problems. The L-shaped method suggested by Van Slyke and Wets was the first such example [92] and applied the Benders’ partitioning procedure [9] to stochastic linear programs. Variants of this approach have formed a significant body of algorithmic research for stochastic programming [78]. Convex generalizations of this problem have recently seen much study through the use of both 1
This work has been submitted to the SIAM Journal of Optimization in February 2008.
1
primal-dual and dual decomposition methods [10, 78]. An alternate approach to tackling large probability supports or continuous distributions is through the use of sampling techniques. Such approaches provide confidence bounds on the optimal value (see [23, 22, 45, 39, 47, 81, 82]) (see section 2.4 for a detailed outline of literature). However, there are a multitude of practical problems where convexity cannot be assumed [31, 11, 68]. Later we provide some motivating examples, drawn from financial planning and Stackelberg equilibrium problem amongst others. Existing stochastic programming approaches, particularly those constructed under convexity assumptions, cannot work on such problem classes. This paper concentrates on precisely such a gap; specifically, we provide a scheme for addressing general stochastic nonlinear programs that satisfies three crucial properties. First, the algorithm is globally convergent; second it possesses rapid local convergence properties and third the algorithm scales well with the size of the probability support requiring that it has a certain degree of parallelizability with respect to scenarios. The major contributions of this paper lie in the design of a novel SQP framework for stochastic nonlinear programs. In keeping with the spirit of SQP methods, we solve the original problem through the solution of a sequence of convex quadratic programs (QPs). A sparse quasi-Newton update is proposed that maintains positive-definiteness and therefore ensures that the resulting stochastic QPs are convex. This allows for the employment of QP duality theory to decompose the stochastic QPs into scenario-based problems. A key difference between existing SQP methods and our framework lies in the techniques for solving the resulting QP approximations, which despite being arbitrarily large, possess a decomposable structure. We propose two methods rooted in the ideas of Benders decomposition for solving such structured QPs: (1) an inexactcut L-shaped method [82, 81] and (2) a quadratic extension of the trust-region method of Linderoth and Wright (referred to as LW, hereafter) [56]. By conducting a parallelizable linesearch on an augmented Lagrangian merit function, we ensure global convergence of the method. It turns out that an extension of the SQP method to a Benders framework is not that straightforward. A disadvantage of Benders decomposition is that estimates of multipliers are not directly obtained from the solvers. We propose a scheme that obtains such estimates through the solution of a single low-dimensional quadratic program . The structure of the stochastic nonlinear program dictates that the resulting 2
QP approximation is a stochastic program with random recourse. Such problems are less well studied and using some classical results, we discuss a characterization for the set of feasible first-stage decisions, under such an assumption on recourse. Finally, the performance of the method is examined from the standpoint of scalability and local behavior on a stochastic nonlinear test problem set. We observe that the effort for solving large stochastic NLPs grows linearly with the size of the probability support. In addition, the framework can address problems whose deterministic equivalents are well over a few hundred thousand variables and constraints. It should be noted that significant research effort has been applied in the direction of developing algorithms for multiperiod stochastic programming but primarily in the linear programming domain. Beginning with the work by Birge [12] on nested decomposition methods, a multitude of approaches has been studied for solving such problems, ranging from cutting-plane schemes to interiorpoint methods, amongst others [13, 79]. Given the massive dimensionality of such problems, Monte-Carlo sampling methods have also been employed in an effort to obtain estimators [46, 84]. We contend that this paper demonstrates to some extent that Benders-SQP approaches could conceivably be extended to the multiperiod nonlinear programming domain. The first step in such a direction would require the development of a multiperiod stochastic QP solver. An obvious choice would be an extension of nested Benders techniques to allow for handling quadratic programs. This however is beyond the scope of this paper, the primary focus being on two-period stochastic nonlinear programming. The remainder of the paper is organized into six sections. In section 2, we provide some applications of nonlinear programming under uncertainty and survey a subset of past approaches. Section 2.5 is unique, dealing with some questions pertaining arising from our random recourse setting. In section 3, we introduce our model for stochastic nonlinear programming and present an outline of the SQP-based framework. The two techniques used for solving the resulting QP approximations are presented in section 4. In section 5, we demonstrate the scalability and efficiency of the algorithm on a test problem set of stochastic NLPs. We conclude in section 6.
3
1.2 Background Problems of nonlinear programming under uncertainty arise out of the modeling of decision-making problems in which some problem data is deemed to be uncertain and decisions need to be made prior this uncertainty reveals itself. We provide three motivating problems in sections 2.1–2.3 while others may be found in [13, 80]. In section 2.4, we provide a summary of preceding research in this area. Finally in section 2.5 we discuss some properties of the underlying QP approximation.
1.2.1 Financial Planning We present a modified version of the financial planning problem discussed in [80]. Suppose we have a budget of W0 and n investment opportunities with random returns in the next time period given by R(ω) = (R 1 (ω), . . . , R n (ω)), where ω is a random variable. The decision-making problem requires the allocation of wealth across n assests, as denoted by x = (x 1 , . . . , x n ), so as to maximize the expected utility of wealth. The utility function is concave, nondecreasing and twice differentiable. The resulting optimization problem requires the optimal allocation so as to maximize expected utility, subject to initial wealth constraints and is given by (FinPlan): " Ã
FinPlan
minimize −IEω U x 1 ,...,x n
subject to
n X
n X
!#
(1 + R i (ω))x i
i =1
x i = W0 .
i =1
1.2.2 Newsvendor Problem Our next example is that of the newsvendor problem, a nonlinear analog of what is commonly seen in stochastic linear programming literature [56, 13]. On a given day, a newsvendor buys x newspapers, at a nonlinear purchase price (this may arise due to, say, discounting by the supplier) c(x) before the demand for newspapers is known. The newspapers are sold after the materialization of demand, resulting in a concave revenue function q(y) for sales y and the unsold newspapers, w, are returned back to the supplier at a rate r . If we assume risk neutrality of the newsvendor, then the objective is to maximize the
4
newsvendor’s profit. Note that having the option of returning the newspapers at a positive price means that maximizing scenario-wise earning is itself an optimization problem with decision variables y and w. The nonlinear newsvendor problem then becomes £ ¤ ω ω minimize xc(x) − IE q(y ) + r w ω ω
NV
x,u ,v
yω + wω = x y ω ≤ dω
subject to
x, y ω , w ω ≥ 0,
∀ω ∈ Ω.
In both of the earlier examples, either the allocation of wealth (FinPlan) or the number of newspapers (NV) are here-and-now decisions; such decisions are required to be made before the uncertainty reveals itself.
1.2.3 Stochastic Equilibrium Problems Our third source problem falls within the class of stochastic equilibrium problems. Such decision-making problems require a first-stage decision x prior to the specification of a recourse decision y ω in accordance with a scenariospecific equilibrium problem parametrized by x. This may be viewed as one variant of a Stackelberg setting [88] in which a leader makes a decision x, after which a set of followers compete in a Nash fashion. In such a setting, the leader would solve the following optimization problem: StochEquil
£ ¤ minimize f (x) + IE h ω (y ω ) ω x,y
y ω ∈ SOL(F ω (x), K ω (x)), subject to
∀ω ∈ Ω,
c(x) = 0, x ≥ 0,
where f and h ω represent the cost making decision x and (recourse) decisions y ω , respectively. Furthermore, SOL(F, K ) represents the set of solutions u to the variational inequality [27] (y − u)T F (u) ≥ 0,
∀ y ∈K,
where K is a closed convex set.
5
Such problems represent stochastic generalizations of mathematical programs with equilibrium constraints or MPECs [59]. While such problems do not satisfy conventional regularity conditions at any feasible point, a common approach has been through the solution of a sequence of regularized (and therefore wellposed) nonlinear programs [75].2 Similar approaches for solving stochastic MPECs would require efficient and scalable stochastic NLP solvers, a focus of this paper. In fact, extensions of SQP methods to address ill-posed NLPs such as MPECs have been suggested by Anitescu [3].
1.2.4 Summary of Literature To understand past approaches, some observations about typical stochastic programs such as the newsvendor problem are worth noting. The problem below is a stochastic linear program. StLP
£ ¤ minimize c T x + min q ωT y ω x
Ax = b subject to
A ω x + B ω y ω = bω x, y ω ≥ 0,
∀ ω ∈ Ω.
Under the assumption of discreteness of ω, the stochastic program can be formulated as a deterministic equivalent program (DEP). Specifically, given that the parameters are known under all realizations of the given distribution, it is ¡ ¢ theoretically possible to solve the DEP directly in x, y 1 , y 2 , ..., y |Ω| . Practically, the size of the probability support, namely |Ω|, can be arbitrarily large making such a solution technique unrealistic. Some early algorithmic efforts have studied schemes that take advantage of the structure of the resulting problems [13]. When the resulting problem is a stochastic linear program, Benders’ decomposition can be employed to split the DEP into a master problem and a set of |Ω| subproblems. The resulting method is the L-shaped method of Van Slyke and Wets [92] (see chapter 5 in [13] for a detailed explanation of such a method). Whenever well defined, we define Q(x; ω) as the optimal value of the problem (subω ), and the recourse function Q(x) = IE[Q(x; ω)]. The master and 2
Regularity, in this context, implies that the resulting NLP satisfies conventional regularity conditions.
6
subproblems are defined as master
subω minimize q ωT y ω ω
minimize c T x + Q(x) x
subject to
Ax = b,
y
and subject to
x ≥ 0,
A ω x + B ω y ω = bω y ω ≥ 0,
respectively. The original L-shaped method assumed B ω to be fixed with respect to ω (called fixed recourse). The problem (StLP) does not have this property and instead has random recourse. The issue of random recouse is discussed more deeply in section 2.5 while more general questions are considered in [54]. The basic idea of the L-shaped method is to approximate the nonlinear term
Q(x) by a set of linear constraints (which bound the recourse function from below), called optimality cuts [13]. It can be shown that for a stochastic linear program, Q(x) is a piecewise linear convex function of x [13, 94]. Optimality cuts manifest themselves as linear inequality constraints in the master problem and together form an outer linear approximation of Q(x). Hence in every iteration of the L-shaped method, the master problem solves for the minimum of the sum of c T x and the outer linear approximation of Q(x), subject to the original constraints. Convergence of a scheme is shown in [13] while an inexact generalization is proved in Zakeri et al. [98] and Au et al. [4]. Although significant research has been carried out on stochastic linear programs [46, 38], stochastic convex programs or more generally stochastic nonlinear programs have been relatively understudied. Perhaps the earliest algorithmic work on stochastic convex programs is by Ziemba [102] which spoke of methods using steepest descent or using the explicit expression of the DEP. Theoretical work in explicitly stating the KKT conditions for the stochastic convex program with ω as a continous random variable is found in the work by Rockafellar and Wets [77]. Stochastic convex programs with linear and bound constraints were solved by an interior point method by Blomvall [14] and Blomvall and Lindberg [15] by assuming a special separable problem structure. Second order primal-dual approaches for solving stochastic mathematical programs with equilibrium constraints (MPECs) [59] were suggested by Shanbhag et al. [81, 64]. Such methods ensure convergence to second-order KKT points. Similar convergence behavior is seen in the barrier-Lagrangian scheme suggested by Zhang and Shanbhag [99].
7
One of the difficulties in constructing algorithms for stochastic nonlinear programs is that the value function of the second-stage decisions is not differentiable with respect to the first-stage variable. The last few years has seen a significant increase in Newton–barrier methods for solving stochastic convex programs [101, 100]. The methods use log–barrier functions to force interior solutions to the second stage problems, making it possible to compute gradients and hessians of the value function. But this comes at the expense of going through another loop of barrier parameter updates. In a recent paper by Liu and Zhao [57] multistage stochastic nonlinear programs were addressed by using an SQP method with the scenario analysis technique of Rockafellar and Wets [76]. Their algorithm reformulates the objective of (SNLP) into a sum of scenario based objectives each with a different first-stage variable. The nonanticipativity condition, which requires that all states which are observationally undistinguishable should have identical optimal value for the decision variable, is imposed explicity by adding a set of linear equality constraints.
1.2.5 Random Recourse It is notable that our problem possesses random recourse. Even for the case of linear recourse, such a recourse model has been seen as an impediment to the construction decomposition algorithms like the L-shaped method. Since this type of recourse is not commonly used in literature, we provide some insight into the issues that arise due to random recourse. We refer the reader to [13], page 87, for a short discussion on why random recourse is avoided. Our main result is an extension of a condition (called the W–condition) first suggested in [94]. This condition specified that under some continuity conditions, one could obtain a clear characterization of the set of first-stage decisions that resulted in bounded and feasible second-stage problems with probability one. We extend precisely such a result to the QP setting. We limit the discussion of this section to stochastic linear and quadratic programs, since the latter are what are our QP approximations. Consider the
8
problem StP
£ ¤ minimize f (x) + min h ω (y ω ) x
Ax = b subject to
A ω x + B ω y ω = bω x, y ω ≥ 0,
∀ ω ∈ Ω.
where f (·) and h ω (·) are either linear or quadratic functions. Define q ω (u ω ) as the dual function of the second stage problem (P2ω (x))
min { h ω (y ω ) | B ω y ω = b ω − A ω x, y ω ≥ 0}
and Uω as the appropriate set of constraints on the dual variables u ω . The object of concern in stochastic programming is the set of permissible first stage decisions x, such that (StP) is feasible and bounded with probability 1. In fact, the L-shaped method (and our variants) produce solutions over a possibly larger set of realizations resulting in a more restrictive set of feasible first-stage decisions. Indeed x = x ∗ is considered “optimal” for (StP) if it is optimal with probability 1. Since the first stage problem is deterministic, it is enough to consider the set of permissible x as the set K defined as K = {x | (P2ω (x)) is feasible and bounded with probability 1}. L-shaped methods use a master-subproblem structure by solving for x first and then solving (P2ω (x)) for all ω ∈ Ω. Our interest is in finding how K is related to the feasibility and boundedness of each (P2ω (x)), ω ∈ Ω. For all x for which (P2ω (x)) is feasible, we define
Q(x; ω) := min{hω (y ω ) | B ω = bω − A ω x, y ω ≥ 0}
and
Q(x) := IE [Q(x; ω)] .
If (P2ω (x)) is infeasible for some x, we define Q(x; ω) = +∞. We can then rewrite K using the recourse function as K = {x | − ∞ < Q(x; ω) < ∞ with probability 1}. For a characterization of K, we will first concern ourselves with feasibility of
9
(P2ω (x)). In accordance with the work in [94], we define three types of sets feasibility set
K 2 := {x | Q(x; ω) < +∞ with probability 1},
strong feasibility set
K 2s := {x | Q(x) < +∞}, K 2 (ω) := {x | Q(x; ω) < +∞}.
elementary feasibility sets
Our endeavor is to find a characterization of K 2 in terms of the sets K 2 (ω). It follows that K 2s ⊆ K 2
and
\ ω∈Ω
K 2 (ω) ⊆ K 2 .
e to be the smallest relatively closed subset of Ω of measure 1. i.e. Define Ω e= Ω
\
S.
S ⊆ Ω, S closed, IP(S) = 1
e is called the support of the probability measure IP and is a measurable set Ω e the set of recourse maif IP is a regular Borel measure. Let B = {B ω : ω ∈ Ω},
trices corresponding to the smallest scenario set of measure 1. An algorithmic procedure that uses master–subproblem structure outlined in section 1.2.4 encounters the elementary feasibility set. A characterization of K 2 in terms of the elementary feasibility sets K 2 (ω) is essential for the development of decomposition algorithms. The following theorem sourced from [94] provides us an important characterization of K 2 . Theorem 1.2.1 (Walkup and Wets, theorem 3.7 [94]) If the restriction of the setvalued map pos(·) to B is continuous, we have e = K 2 = {x | Q(x; ω) < +∞ ∀ ω ∈ Ω}
\
K 2 (ω).
e ω∈Ω
We denote the restriction of a set-valued map F : X 7→ Y to a set K ⊆ X by F |K where X , Y are metric spaces. Specifically F |K is defined as [5]
F |K =
F (x)
if x ∈ K
;
if x ∉ K .
For a finite set Ω, theorem 1.2.1 implies that K 2s = K 2 . This still does not elim˜ Q(x) ˜ < ∞ but Q(x; ˜ ω) e = +∞ for some inate the possibility that for some x = x, 10
e ∈ Ω. i.e. we have ω
K2 =
\
K 2 (ω) ⊇
e ω∈Ω
\ ω∈Ω
K 2 (ω).
e 6= Ω, i.e. IP(ω) e = 0 for ω e as above. When Ω is finite such a situation arises when Ω
Essentially, such scenarios would introduce constraints that further restrict the feasible region. While on paper such an x˜ would qualify as feasible for the stochastic program, it could create numerical discrepancies in the evaluation ˜ ω)]. It is left to prue × Q(xe; ω) e = 0 × ∞ as part of the summation IE[Q(x; of IP(ω) dence of the programmer whether he chooses to solve the stochastic program e or find his way around possible numerical difficulties. for the restricted set Ω This is also an important insight into decision making under uncertainty. Often in real life problems, like (FinPlan) for instance, scenarios that lead to infinite returns also have an infinitesimal likelihood of occurence. What we have presented above is a mathematical justification of why decision making uncertainty should be done on the smallest scenario set of measure 1. In fact, the L-shaped method may result in an no feasible first-stage solution if one added appropriately chosen constraints corresponding to measure zero events. Next, we discuss the boundedness of second stage problems. Assessing the boundedness of the second stage problems is equivalent to assessing the feasibility of their duals. Define analogously three dual feasibility sets K 2∗ := {x | Q(x; ω) > −∞ with probability 1}
dual feasibility set
K 2s∗ := {x | Q(x) > −∞}
dual strong feasibility set elementary dual feasibility sets
K 2∗ (ω) := {x | Q(x; ω) > −∞}.
In the light of these definitions, the set of permissible x for (StP), K is K = K 2 ∩ K 2∗ .
(1.1)
We now desire a characterization of K 2∗ in terms of K 2∗ (ω). Since the duals of the second stage problems are involved here, this characterization depends on whether h ω (·) is linear or quadratic. For the linear case we have the following result from [94]. Theorem 1.2.2 (Walkup and Wets, corollary 3.9 [94]) For a stochastic linear program, if the restriction of the set-valued map pos(·) to {[B T , −B T , I ] : B ∈ B } is
11
continuous we have e = K 2∗ = {x | Q(x; ω) > −∞ ∀ ω ∈ Ω}
\
K 2∗ (ω).
e ω∈Ω
Theorems 1.2.1 and 1.2.2 together constitute the W–condition for stochastic linear programs [94]. Definition 1.2.3 (W–condition for stochastic linear programs) A stochastic linear program with a set of recourse matrices {B ω : ω ∈ Ω} is said to satisfy the W– condition if the restrictions of set-valued maps pos(B ω ) and pos([B ωT , −B ωT , I ]) to B and {[B T , −B T , I ] : B ∈ B } respectively are continuous. In [94], sufficient conditions are provided for stochastic linear programs to satisfy the W–condition. Theorem 1.2.4 (Walkup and Wets theorem 3.7 [94]) Any one of the following is a sufficient condition for the stochastic linear program to satisfy the W–condition e pos(B ω ) is a pointed cone and no column of B ω has zero 1. For each ω ∈ Ω,
norm. e pos(B ω ) is a (not necessar2. There exists an integer k such that for all ω ∈ Ω,
ily fixed) subspace of ℜm of dimension k. 3. B ω is constant throughout Ω. The third alternative in the above result is the familiar fixed recourse condition seen in a variety of stochastic programming literature and its origin lies in the satisfaction of the W–condition. We now turn to the quadratic case, where we restrict ourselves to strictly convex functions h ω (·). The contribution of this section is in indentifying an analogous W–condition for such stochastic quadratic programs. An important observation is that a feasible strictly convex quadratic program is always bounded below (by its unconstrained minimum). This is a point of departure from LPs where feasibility does not necessarily imply boundedness. i.e. for strictly convex QPs we have K 2 (ω) ⊆ K 2∗ (ω).
12
So if the condition in theorem 1.2.1 is met, K2 ⊆
\
K 2∗ (ω) ⊆ K 2∗ .
e ω∈Ω
Using theorem 1.2.1 we get the following theorem. Theorem 1.2.5 For a strictly convex stochastic quadratic program, if the restriction of pos to B is continuous, then K = K2. We conclude that for such programs, it is sufficient to provide a characterization of K 2 in terms of K 2 (ω). In other words, theorem 1.2.1 alone constitutes a W–condition for strictly convex stochastic quadratic programs. Definition 1.2.6 (W–condition for strictly convex stochastic quadratic programs) A strictly convex stochastic quadratic program with a set of recourse matrices {B ω : ω ∈ Ω} is said to satisfy the W–condition if the restriction of pos(·) to B is continuous. The following theorem from [94] provides a sufficient condition for the restriction of the pos mapping to be continuous. Denote L (C ) to be the lineality space of a cone C defined by
L (C ) = C ∩ −C . Theorem 1.2.7 (Walkup and Wets, theorem 2 in [93]) Suppose Z is a subset of ℜp×n , k is an integer, and for every matrix Ab ∈ Z , b = k, (a.) dim L (pos( A))
(b.) there exists a neighborhood N of Ab such that if any column Abi of Ab lies in b then the corresponding column A i of any matrix A ∈ N ∩ Z L (pos( A)), lies in L (pos(A)). Then the restriction of pos to Z is continuous. The set B forms a subset of ℜp×n . It can be verified that theorem 1.2.4 is an application of theorem 1.2.7. An extreme case of theorem 1.2.7 is when for any
13
b is a pointed cone and none of the columns of Ab are zero. This imAb ∈ Z , pos( A) b = 0, making condition (b.) of theorem 1.2.7 vacuously true plies dim L (pos( A))
and the pos mapping continuous on Z . An even more extreme case is when Z is a singleton (fixed recourse) which makes both (a.) and (b.) vacuously true. In this light, conditions of theorem 1.2.4, including the popular fixed recourse can be seen as too strong for the W–condition to be satisfied. We now show that for a finite set Ω, condition in definition 1.2.6 can ensured simply. Before proving the continuity of the set-valued map, we recall the definition of upper semicontinuity in the context of such maps. Definition 1.2.8 (Upper semicontinuity of set-valued maps [5]) Let X and Y be metric spaces. A set–valued map F : X 7→ Y is called upper semicontinuous at x ∈ dom(F ) if and only if for any neighborhood U of F (x), ∃ η > 0 such that ∀ x 0 ∈ B X (x, η),
F (x 0 ) ⊂ U,
where B X (x, η) is a ball defined on X centered at x and with radius η. A map is said to be upper semicontinuous if it is upper semicontinuous for all x ∈ X . Proposition 1.2.9 For a finite set Ω, the restriction of pos(·) to B is continuous. Proof : Proving such a result requires proving the upper and lower semicontinuity of the map. In general the pos map when considered as a (cone–valued) mapping from ℜp×n to ℜp is lower semicontinous [93]. Hence for definition 1.2.6 to hold, we need only to prove the upper semicontinuity of pos restricted to B . Denote the restriction of pos to B by pos |B . Let B ω ∈ dom(pos |B ) e and and consider a neighborhood U around pos(B ω ). Since Ω is a finite set, Ω
B are also finite. By the finiteness of B , there exists a neighborhood N ⊂ ℜp×n of B ω such that N ∩ B = {B ω }. Furthermore, for all B ∈ N\B ω , pos |B (B ) = ;. So for all B ∈ N, pos(B ) ⊂ U. We now summarize the contents of this section. The W–condition is needed to get a characterization of K 2 in terms of K 2 (ω) and K 2∗ in terms of K 2∗ (ω). We proved that for strictly convex stochastic quadratic programs it is enough to get
14
a characterization of K 2 in terms of K 2 (ω). If the restriction of pos() to B is continuous, we get K2 =
\
K 2 (ω).
e ω∈Ω
So for strictly convex stochastic QPs we say, that if the restriction of pos to B is continuous, the W–condition is satisfied. But for a finite set Ω, the restriction is continuous implying that the W–condition for stochastic QPs holds automatically. By observing definition 1.2.3 we see that for finite Ω, the W–condition for stochastic LPs also holds automatically. Finally, no extra assumptions on recourse matrices are needed for stochastic LPs and strictly convex stochastic QPs when one assumes that Ω is a finite set. A final benefit of this analysis is that e ensuring that the it leads us to assumption 1.3.1(e) which ensures that Ω ≡ Ω, L-shaped method works with appropriate set of feasible first-stage decisions.
1.3 An SQP Framework Our model for stochastic nonlinear programming falls under the category of recourse programs, and is specifically a two-stage recourse program. Such programs are characterized by two decisions [13]: 1. a first-stage decision taken before the uncertainty is resolved 2. a recourse decision taken after the uncertainty reveals itself to help compensate for the misjudgement made in the first-stage decisions. Some stochastic programming models capture the dependence of the recourse decision on the first-stage decision by explicitly including the first stage variable as an argument of the second-stage objective functions [57]. In our model, we keep the the objective function dependent only on the second-stage variables and capture this dependence through the second-stage constraints, which we allow to be nonlinear. As we see, this provides the problem with a nearly separable structure that proves useful in algorithm design (see [13] or Ruszcynski and Shapiro [80] for a review of stochastic programming models). Various recourse models exist in literature defined for linear constraints. We refer the reader to ongoing work by the authors [54] for an account of recourse in nonlinear and competitive settings.
15
1.3.1 Sequential Quadratic Programming We make the following assumptions about (SNLP): Assumption 1.3.1
(a.) The optimal value of (SNLP) is bounded and denoted
∗
by v . (b.) The functions f (·), h ω (·), u(·), a ω (·) and d ω (·) are nonlinear continuously differentiable multifunctions for all ω ∈ Ω. (c.) For all ω ∈ Ω the Jacobian of d ω (·), denoted J dω (y ω ) is full row-rank. (d.) For every x satisfying u(x) = 0 and x ≥ 0, and ∀ ω ∈ Ω, there exists y ω satisfying a ω (x) + d ω (y ω ) = 0 and y ω ≥ 0. (e.) For all ω ∈ Ω, IP(ω) > 0. Assumptions (a), (b), (c) are standard assumptions made for any SQP method. (d) is the nonlinear analog of relatively complete recourse [13]. Assumption (e) implies that Ω is the smallest relatively closed subset of Ω with measure 1. In e Note that (SNLP) requires the satisfaction of constraints arising effect, Ω ≡ Ω. from each ω ∈ Ω. However the objective of (SNLP) contains only those functions corresponding to scenarios for which IP(ω) > 0. Without (e) any stochastic program could be appended with any number of dummy constraints without a change in the objective value by posing them as having originated from scenarios of zero probability. Assumption (e) prevents such possibilities and keeps the problem well defined. The SQP approach to solve nonlinear optimization problems found inception in the Ph.D. thesis of Wilson in 1963 [96]. But the method became truly popular after Hessian-update methods for problems with indefinite Lagrangian Hessians were studied by Mangasarian and Han [16, 34, 35]. Through several papers by Powell in the late 1970s [74, 72, 73], Han’s theorems on the convergence of SQP methods became widely known in the optimization community. The SQP approach is an iterative method with the basic idea being to model the nonlinear program at a given iterate x k by an approximate quadratic program. The solution of this QP provides a direction p k while an acceptable stepsize αk is obtained by conducting a linesearch on an appropriately defined merit function. The new iterate x k+1 may be defined as x k + αk p k . Under some assumptions, it can be shown that the sequence of iterates converges to the 16
KKT point of the original nonlinear program. SQP methods offer great flexibility in the construction of the QP approximations in order to be able to deal with regions with an indefinite Hessian, with the overarching goal of driving the sequence of solutions to a KKT point. With an appropriate choice of Hessian matrices, SQP methods can be shown to be equivalent to a Newton or quasiNewton method applied to solving the KKT system of the original nonlinear problem [68]. As a result of this equivalence, in regions close to the KKT point of the orginal problem, with exact Hessians, SQP methods show similar behaviour to Newton methods – that of quadratic local convergence. If the SQP method is augmented with a merit function, a reduction in which directs the sequence x k toward a KKT point, the SQP method can be proved to be globally convergent [16]. There are several practical SQP methods tailored for different kinds of applications; NPSOL [30] for dense problems, SNOPT [29] for sparse problems and
filterSQP – a filter method [28] – being some of them. The method that follows is the one employed in the solver SNOPT by Gill, Murray and Saunders [29], which we shall apply to (SNLP). We now present a quick overview of SNOPT.
1.3.2 A Linesearch SQP Method We provide a description of our SQP method which is based on the linesearch SQP method in [29]. Consider a general nonlinear program GNP
minimize f (x) x
subject to
c(x) = 0,
:λ
x ≥ 0.
:µ
If (x∗ , λ∗ , µ∗ ) is a local minimum of (GNP) and s∗ is the slack variable corresponding to nonnegativity bounds on x, it satisfies the first-order KKT conditions given as ∇ f (x∗ ) − (λ∗ )T ∇c(x∗ ) − µ∗ = 0 c(x∗ ) = 0 0 ≤ s∗ ⊥ µ∗ ≥ 0 x∗ − s∗ = 0.
17
We define the linearization of the constraint at x k as c L (x; x k ) := c k + ∇c k (x − x k ), and the departure from linearity at x k as d L (x; x k ) := c(x) − c L (x; x k ). The modified Lagrangian L is defined as
L(x; x k , λk ) := f (x) − λTk dL (x; x k ). Observe that L(x k ; x k , λk ) = f k and ∇L(x k ; x k , λk ) = ∇ f (x k ) := g k . This SQP method solves (GNP) as sequence of QPs, each of which minimizes a quadratic approximation of the modified Lagrangian with respect to the linearization of the constraints. So at (x k , λk ), the QP solved is QPk
minimize f (x k ) + g kT (x − x k ) + 12 (x − x k )T Hk (x − x k ) x
subject to
c k + ∇c k (x − x k ) = 0 x ≥0
: λ¯ k : µ¯ k ,
where Hk is a positive definite approximation to ∇2 L(x k ; x k , λk ). Let (x¯k , λ¯ k , µ¯ k ) be a solution to (QPk ) and s¯k be the slack variable. The optimality conditions for (QPk ) are g k + Hk (x¯k − x k ) − λ¯ Tk ∇c k − µ¯ k = 0 c k + ∇c k (x¯k − x k ) = 0 0 ≤ x¯k ⊥ µ¯ k ≥ 0 x¯k − s¯k = 0. It can be seen that the current solution (x¯k , λ¯ k , µ¯ k , s¯k ) can be regarded as an approximation to (x∗ , λ∗ , µ∗ , s∗ ). Algorithms for nonlinear optimization problems are faced with the twin goals of improving either feasibility and optimality at each iterate. The merit function combines both these concerns within a single function. In moving from one iterate to another, the new iterate is determined by conducting a linesearch on the augmented Lagrangian merit function, which is defined as
M(x, λ, µ, s) = f (x) − λT c(x) − µT (x − s) + 12 ρkc(x)k2 + 21 ρkx − sk2 , where ρ is the penalty parameter. The new step is searched for along the direc-
18
tion joining of the current iterate to the new solution
x¯k − x k λ λ λ¯ − λ k k+1 k k = + αk . µk+1 µk µ¯ k − µk s k+1 sk s¯k − s k
x k+1
xk
Descent in the merit function is needed to ensure globally convergent behavior. The penalty parameter is raised in the case when sufficient descent is not obtained. In fact, we prove that given a solution to (QPk ), there exists a finite ρ for which descent with respect to the merit function can be guaranteed. We end this subsection with an introduction of the SQP framework to the case of stochastic NLPs: SNLP
£ ¤ minimize f (x) + IEω h ω (y ω ) x,y
u(x) = 0 subject to
a ω (x) + d ω (y ω ) = 0 x, y ω ≥ 0
: λx : p ω λω : µx , p ω µω ,
∀ω∈Ω
where λx , p ω λω , µx and p ω µω denote the Lagrange multipliers for the respective constraints scaled by the scenario probability mass for convenience. Following are some notational rules we follow during our description of the SQP algorithm. Some of these notations might change meaning later in the paper; the new meaning will be made clear as and when necessary. • We denote y, z, λ, µ, ν and c(z) as
y :=
y1
λx
µx
u(x)
a (x) + d (y 1 ) p µ1 p λ1 Ã ! Ã ! 1 1 1 1 y2 λ x 2 2 2 a 2 (x) + d 2 (y ) p2µ p2λ , ν := and c(z) := , µ := , z := , λ := .. . µ y . . . . .. .. .. y |Ω| |Ω| |Ω| |Ω| a |Ω| (x) + d |Ω| (y ) pωµ pωλ
• Subscripts and superscripts: – As with the Lagrange multipliers of (SNLP) above, we use a superscript x to denote terms associated with the cost function f (x) or 19
the constraint u(x). second-stage variables, constraints or objectives will be denoted by ω ∈ {1, . . . , |Ω|} as superscript. All quantities that are given as a part of the definition of SNLP – viz., a ω (.), h ω (.) etc. use the scenario number as subscript and thus are distinguished from quantities that we introduce into the problem – viz. y ω , λω etc. – As an exception to the above rule, the expectation symbol when subscripted with ω as “IEω ” means “expectation with respect to random variable ω”. – We reserve the subscript k to denote the SQP iteration to which the subscripted variables belong: e.g. x k , y kω , λkx etc. For a function F (x), F k will denote its value at x k , unless stated otherwise. – The superscripted asterisk denotes the optimal for (SNLP): (z ∗ , ν∗ ) is the KKT point of (SNLP). – We use the subscript L and Q along with the function to depict its linearization and its quadratic approximation respectively. • Jacobians and Hessians (or approximations to the Hessian, whichever is applicable) will be denoted by J and H respectively. They will be affixed with suitable subscripts and superscripts to denote which SQP iteration and functions they pertain to. • Let δL denote the departure from linearity of c(z) at the current iterate z k while c L denote its linearization. δL (z; z k ) := c(z) − c L (z; z k ), where c L (z; z k ) := c(z k ) + J kc ∆z k . • Lastly, for any variable x, we use ∆x k to denote the difference (x − x k ).
1.3.3 Constructing QP Approximations SQP literature commonly refers to such approximations as QP subproblems [29]. We will continue to refer to them as approximations and instead reserve the term subproblem for scenario-wise problems obtained from decomposing the stochastic program as stated in section 1.2.4. The modified Lagrangian is
20
defined as
L(z; νk , z k ) := f (x) + IEω hω (y ω ) − λTk δL (z; z k ). £
¤
£ ¤ b (z; νk , z k ) := f (x) + IEω h ω (y ω ) − We distinguish it from the true Lagrangian L
λTk c(z) − µTk z. Let g k := ∇Lk . By the procedure described in the preceding section, the QP approximation is StQPk
minimize Lk + g kT ∆z k + 21 ∆z kT Hk ∆z k z
subject to
c L (z; z k ) = 0 z ≥ 0,
where Hk  0 is a positive definite approximation to ∇2 Lk . Continuing the same notation as for (SNLP), let (z¯k , ν¯k ) be the KKT point of (StQPk ). We assume the following about the nature of the solution of (StQPk ). Assumption 1.3.2 At the solution of (StQPk ), (z¯k , ν¯k ), the constraint Jacobian J kc is full row rank and the Hessian of the Langrangian ∇2 L(z¯k , ν¯k ) is positive definite in the tangent space of the constraints. The modified Lagrangian function L may be written as L = Lx + IEω [Lω ] , where
Lx (x; x k , νk ) = f (x) − λkxT u L (x) − IEω λωT k a ωL (x; x k ) £
¤
ω Lω (y ω ; y kω , λω k ) = h ω (y ).
From a structural standpoint, the Hk needs to satisfy the following conditions: 1. It must have a block diagonal structure with blocks corresponding to x and each y ω , ω ∈ Ω and 2. each such block should be positive definite. We will return to the question of constructing such Hessians, with some added requirements in section 1.3.5. Given Hessians of such structure, the quadratic approximations of the Lagrangian functions may then be stated as
LQx (x; x k , νk ) = Lkx + ∆x kT ∇Lkx + 12 ∆x kT Hkx ∆x k ω ω T ω ω T ω ω 1 LQω (y ω ; y kω , λω k ) = Lk + ∆(y k ) ∇Lk + 2 ∆(y k ) Hk ∆y k .
21
With such a structure, the QP subproblem (StQPk ) can now be written as h i x ω ω ω ω StQPk minimize LQ (x; x k , νk ) + IEω LQ (y ; y k , λk ) x,y
subject to
J ku ∆x k = −u k
: λ¯ x
J k ω ∆y kω + J k ω ∆x k = −a ω,k − d ω (y kω ) : p ω λ¯ ω d
a
x, y ω ≥ 0,
: µ¯ x , p ω µ¯ ω
∀ω ∈ Ω.
In section 4, we provide two approaches for solving such a class of problems: • The inexact L-shaped method • A trust region based L-shaped method, which is a QP extension of the trust-region method of LW [56]. Both are based on using Benders’ decomposition for quadratic programs. The inexact L-shaped method uses inexact optimality cuts to model the recourse function. The inexactness is sequentially reduced to make cuts near-exact around the optimal x. Moreover, the inexactness allows for the solution of dual problems to feasibility (as opposed to optimality) alleviating the computational burden significantly. The trust region method uses exact cuts but confines the cuts added to a box-shaped region, which facilitates a close modeling of the recourse function. In the following sections we will assume for the moment that (StQPk ) is solvable and that its KKT point has been obtained.
1.3.4 A Parallelizable Linesearch The augmented Lagrangian merit function for an SQP method on (SNLP) is
M(z, ν, s; ρ) = f (x) + IEω [hω (y ω )] − λT c(z) − µT (z − s) + 12 ρ
X
kc i (z)k2 + 21 ρkz − sk2 .(1.2)
The linesearch procedure is as follows. Having obtained the KKT point (z¯k , ν¯k ) of (StQPk ), we find the slack variables s¯k corresponding to z¯k . The search direction (p k , ξk , q k ) in the primal variables, the multipliers and slacks is constructed as ((z¯k − z k ), (ν¯k − νk ), (s¯k − s k )). For a given ρ, αk denotes the steplength satisfying the Armijo conditions for M(z k + αk p k , νk + αk ξk , s k + αk q k ; ρ). The fol-
22
lowing condition ensures sufficient descent [32]:
pk
T
ξ ∇Mk < − 1 p T Hk p k ≤ 0. k 2 k qk
(1.3)
If (1.3) is not met, we increase the parameter ρ so as to satisfy it. The following lemma shows that this can indeed be done. Lemma 1.3.3 At a given iterate (z k , νk ) there exists a finite ρ¯ such that for all ¯ the sufficient descent condition (1.3) is satisfied. ρ ≥ ρ, Proof : The derivative of M from Eq (1.2) at (z k , νk , s k ) is given by
∇Mk
=
g k − (J kc )T λk − µk + ρ(J kc )T c k + ρ(z k − s k ) −c k −(z k − s k ) µ − ρ(z k − s k )
.
Therefore, we have that
pk
T
à ! c k ξ ∇Mk = g T p k + p T J cT (ρc k − λk ) − ξT k k k k k zk − sk qk
(1.4)
+(q − p)T (µk − ρ(z k − s k )). Observe that since s¯k = z¯k , (q k − p k ) = (z k − s k ). And since z¯k is feasible for (StQPk ), J kc p k = −c k . So using Eq (1.4), condition (1.3) becomes − 21 p kT H p k ρ > Case 1:
− g kT p k − c kT λk +
à !T Ã λ¯ k − 2λk ck
µ¯ k − 2µk
kz k − s k k2 + kc k k2
zk − sk
!
.
(1.5)
The current iterate is infeasible or kz k −s k k2 +kc k k2 > 0. Since (z k , νk ) 6=
(z ∗ , ν∗ ), the right hand side of (1.5) is well defined. Hence there exists a ρ < ∞ such that condition (1.3) is satisfied. Case 2:
The current iterate is feasible or kz k −s k k2 +kc k k2 = 0. Therefore, from
23
(1.4), we have that
pk
T
à ! c k ξ ∇Mk = g T p k + p T J cT (ρc k − λk ) − ξT + (q k − p k )T (µk − ρ(z k − s k )) k k k k k zk − sk qk
= g kT p k + p kT J kcT (−λk ) + (q k − p k )T (µk ). But since c k = 0, we have that J kc p k = 0. Also (q k − p k ) = (z k − s k ) = 0, leaving us to prove that g kT p k < − 12 p kT Hk p k . To prove this, observe that z = z k is feasible for (StQPk ), and at the solution of (StQPk ), ∆z k = p k . Denote the objective of (StQPk ) by F (z). Hence F (z k + p k ) < F (z k ), implying precisely that g kT p k < − 21 p kT Hk p k . Given an αk , the (k + 1)t h iterate is computed as Ã
z k+1
νk+1
!
Ã
=
zk
νk
!
Ã
+ αk
z¯k − z k
!
ν¯k − νk
.
The merit function can be written as
M(z, ν, s; ρ) := f (x) − (λx )T u(x) − (µx )T (x − s x ) + 21 ρ
X i
ku i (x)k2 + 12 ρkx − s x k2
£ ¤ + IEω h ω (y ω ) − (λω )T (a ω (x) + d ω (y ω )) − (µω )T (y ω − s ω ) " # X + IEω 12 ρ ka ω,i (x) + d ω,i (y ω )k2 + 12 ρky ω − s ω k2 . i
The merit function evaluation can be parallelized to a large extent. Specifically, it can be decomposed into first-stage and scenario-wise second-stage evaluations.
1.3.5 A Sparse Quasi-Newton Update The entire analysis that ensued in the previous section relied on the separability of the objective of (StQPk ). Indeed, our solution method for (StQPk ) is also incumbent on this property.
24
The complete requirements from the Hessian approximations are (1.) a block diagonal structure, (2.) block-wise positive semidefinite, with the goal of ensuring local superlinear convergence and allow updating in parallel. A plethora of quasi-Newton updates exist in literature that equip the SQP approach with the convergence properties we desire. In the SQP method used in SNOPT [29], the authors have employed a modified BFGS update . However if we choose to apply this directly on Hk (the approximation to ∇2z L) it does not leave Hk+1 with the block diagonal structure that we need. Sparsity preserving quasi-Newton updates have been the subject of research ever the since the idea of updates originated. Shanno has a sparse variety of the BFGS update [83], there is a sparse PSB update given by Toint [91] and also an explicit quasiNewton update by Lucia [58]. However it is not obvious from the structure of these updates if they can be constructed in a parallelizable manner. Since one of our focii is creating an algorithm suited to parallelization, we do not employ any of the conventional sparsity preserving updates. Quasi-Newton updates can be often be thought of as solutions to optimization problems. In [31], Gill et al. have mentioned this equivalence for the Powell-symmetric-Broyden update. Nocedal et al. [68] have shown the BFGS update to be the solution of the (BFGS)
minimize kW 1/2 (B − B k )W 1/2 kF B
subject to
B pk = rk
(secant condition)
B = BT ,
where p k = z k+1 −z k . If we proceed along the direction of SNOPT [29], we specify r k as b (z k+1 ; z k , νk+1 )−∇L b (z k ; z k , νk+1 ). r k = ∇L(z k+1 ; z k , νk+1 )− ∇L(z k ; z k , νk+1 ) = ∇L
It can be seen from this problem that the objective is not separable – a testimony to the fact that the BFGS update does not retain the block diagonal structure of Hk . The norm in objective is called the weighted Frobenius norm and the weight matrix W is any matrix that satisfies W r k = p k [68]. The unique so-
25
lution to (BFGS) is B k+1 = B k +
r k r kT r kT p k
−
(B k p k )(B k p k )T p kT B k p k
.
(1.6)
Superlinear convergence in a quasi-Newton method for constrained optimization requires the update B k to satisfy the following test [17, 89, 68]: Theorem 1.3.4 ([17]) Suppose assumption 1.3.2 holds and the iterates z k obtained from an SQP method to solve (SNLP) converge to z∗ . Then z → z∗ superk
linearly if and only if lim
b ∗ )(z k+1 − z k )k kP k (B k − ∇2z L
k→∞
kz k+1 − z k k
= 0,
(1.7)
£ ¤−1 c where P k = I − (J kc )T J kc J kcT J k is the projection on the tangent space of the b denotes the true Langrangian (not the modified Lagrangian). constraints and L
Theorem 1.3.4 is an extension to constrained optimization of the original Dennis Moré characterization theorem for quasi-Newton methods [24, 68]. Theorem 1.3.5 (Dennis-Moré) Suppose the iterates z k are generated using a rule z k+1 = z k − B k−1 ∇Fk , for some function F(z) and that z k converges to z∗ . Then z k → z∗ superlinearly if and only if k(B k − ∇2 Fk )(z k+1 − z k )k = 0. k→∞ kz k+1 − z k k lim
(1.8)
It is known that if B k in Eq (1.8) is the sequence of BFGS updates applied in a quasi-Newton method to minimize F(z), then the sequence of iterates {z k } converge superlinearly to z ∗ . The reader may see Nocedal and Wright [68] for a proof. We now define our update rule. Definition 1.3.6 Suppose the Hessian of the Lagrangian function at the kth iterate is denoted by Hk . Then the updated Hessian, denoted by Hk+1 , is given by a relationship, referred to as R(Hk ), where Hk+1 is defined as Hk+1 = R(Hk ) and
26
R(k) is given by x Hk+1 := BFGS(Hkx ), ω Hk+1 := BFGS(Hkω ),
∀ ω ∈ Ω, H0 Â 0,
where BFGS(.) is the BFGS update obtained from problem (BFGS) in Eq (1.6). As opposed to BFGS(Hk ), R(Hk ) preserves the block diagonal structure of Hk . But it is not guaranteed to retain the sparsity structure inside individual blocks. It is a well known property of BFGS updates that they retain positive definiteness as long as p kT r k > 0. We refer the reader to [68] for a proof of this fact. If p kT r k is not sufficiently positive, we choose to simply skip the update [68]. Our results show that this technique works well for our problem. Hence H0  0 =⇒ Hk  0 ∀ k > 0. It is also notable that since the secant condition in problem (BFGS) is separable,
R(Hk ) is feasible for problem (BFGS), but not optimal. We establish now that it nevertheless ensures superlinear convergence. Proposition 1.3.7 Suppose the iterates z k are obtained from the SQP algorithm using update rule R, and suppose they converge to z∗ . Then the update rule satisfies the Dennis-Moré characterisation. Proof : Observe that the KKT conditions for (StQPk ) are Hk p k + g k − J kcT λk − µ = 0 =⇒
b, p k = −Hk−1 ∇z L
b (z, ν) = f (x)+IEω h ω (y ω )−λT c(z)−µT z. Since R is block-wise separable, since L b k and y ω = y ω − (H ω )−1 ∇ y ω L bk, x k+1 = x k − (Hkx )−1 ∇x L k+1 k k
∀ ω ∈ Ω.
(1.9)
b. Consider Eq (1.8) by taking B k = Hk and F = L b k )(z k+1 − z k )k k(Hk − ∇2 L k→∞ kz k+1 − z k k " Ã ! # 2bx x b ω )(y ω − y ω )k X k(Hkω − ∇2 L k(H − ∇ L )(x − x )k k+1 k k k k k+1 k = lim + . k→∞ ω∈Ω kp k k kp k k
lim
27
Take a typical term from the right side like
lim
b ω )(y ω − y ω )k k(Hkω − ∇2 L k k+1 k
kp k k
k→∞
The sequence of iterates {y kω } obtained by the SQP method and the Hessian Hkω b ω . Using the aforementioned property is a BFGS approximation to ∇2 Lω = ∇2 L k
k
of BFGS updates, we find that each such term in the right hand side of this equation is 0. A similar result holds for the term containing Hkx . As a consequence, the result follows. Theorem 1.3.8 Suppose the Hessian updates obtained by rule R are used in our SQP method, and the iterates z k converge to z ∗ . Then z k → z ∗ superlinearly. Proof : Now consider Eq (1.7) with Hk obtained by rule R. 0 ≤ lim
k→∞
b ∗ )(z k+1 − z k )k kP k (Hk − ∇2z L
kz k+1 − z k k
≤ lim kP k k · lim k→∞
k→∞
b ∗ )(z k+1 − z k )k k(Hk − ∇2z L
kz k+1 − z k k
Since {z k } converges, limk→∞ kP k k exists and is finite. By proposition 1.3.7, the rule Hk+1 = R(Hk ) satisfies the Dennis Moré characterization Eq(1.8). This implies that rule R satisfies the condition Eq (1.7) from [17] and z k → z ∗ superlinearly. Importantly, the update R can be applied in parallel on its individual blocks without any information from the other blocks. The termination criteria for SNOPT require that the KKT conditions of (GNP) be satisfied within some accuracy. For prespecified small constants τP and τD define τx = τP (1 + kxk) and τλ = τD (1 + k[λ, µ]k). A point (x, λ, µ) is taken to be a KKT point of (GNP) if it satisfies |c i (x)| ≤ τx , µi ≥ −τλ , |c i (x)λi | ≤ τλ , x i µi ≤ τλ , |[∇ f (x) − λT ∇c(x) − µ]i | ≤ τλ . (1.10) We are now in a position to state the Benders-SQP method (algorithm 1) for solving stochastic nonlinear programs.
28
.
Algorithm 1: Benders-SQP method for Stochastic NLPs 0 initialization k = 1; choose constants τP > 0 and τD > 0, initial first stage decision x 0 , second x ω stage decisions y 0ω , multipliers λ0x , λω 0 , µ0 , µ0 ∀ ω ∈ Ω and initial (SNLP)
Hessian H0 Â 0; while conditions (1.10) are not satisfied do 1
Construct QP approximation (StQPk ) at k t h iterate;
2
Solve (StQPk ) using algorithm ILS or TR;
3
Recover first-stage multipliers by solving problem (MULT) in section 1.4.3;
4
Perform linesearch as shown in section 1.3.4 and get z k+1 ;
5
Construct new Hessian approximation Hk+1 = R(Hk ); end
1.4
Dual Decomposition Methods for the QP Approximations
This section focuses primarily on the decomposition schemes for solving the QP approximations that appear within the SQP framework. After providing some preliminaries on the basic algorithmic framework, namely dual decomposition, we discuss two variants of this decomposition scheme. The first approach uses inexact cuts to approximate the recourse function while the other employs a trust-region to regulate the size of the step made in the master problem. In either case, we consider the following stochastic quadratic program which is in the form of (StQPk ): StQP
minimize x,y
1 T 1 ω T T ω ω 2 x P x + c x + IE 2 (y ) D ω y + d ω y
£
¤
Ax = b subject to
A ω x + B ω y ω = bω x, y ω ≥ 0,
∀ω ∈ Ω.
The matrices P and D ω Â 0 ∀ ω ∈ Ω while B ω is full row-rank ∀ ω ∈ Ω. In addition IP(ω) > 0 ∀ω ∈ Ω. Using the decomposition ideas of Van Slyke and Wets for stochastic linear programs [92], we separate the stochastic QP into a master
29
problem and a set of scenario-based subproblems. Following the notation of section 1.2.4, the master problem (master) and the scenario-specific subproblem (subω ) are defined as minimize x,θ
1 T x P x + cT x + θ 2
y
Ax = b subject to
1 ω T (y ) D ω y ω + d ωT y ω 2
minimize ω and
θ ≥ Q(x)
B ω y ω = bω − A ω x
subject to
y ω ≥ 0,
x ≥0
: vω : uω
respectively. For a positive definite D ω , the Lagrange dual of (subω ) is given by Ã
maximize − 12 ω ω
dsubω
u ,v
uω
!T
Ã
Qω
vω
uω
!
vω
Ã
+ q ωT
uω
!
vω
+ rω
u ω ≥ 0,
subject to where Ã
Qω =
−1 Dω
−1 T Dω B
−1 B Dω
−1 T B Dω B
!
Ã
qω =
,
−1 Dω dω
−1 B Dω d
!
− (A ω x − b ω )
−1 dω . and r ω = − 12 d ωT D ω
If (subω ) is feasible then strong duality holds and the optimal values of (subω ) and (dsubω ) are equal. This allows us to define Q(x; ω) (the random recourse function) in terms of the optimal point of (dsubω ), (u ω∗ , v ω∗ ): Definition 1.4.1 The random recourse function is defined as à !T ω∗ u 1
Q(x; ω) := − 2
v ω∗
Ã
Qω
u ω∗ v ω∗
!
à ! ω∗ u T
+ qω
v ω∗
+ rω
∀ x ∈ K 2 (ω).
The recourse function is defined as the expectation of the random recourse function or
Q(x) := IE[Q(x; ω)]
∀ x ∈ K2.
Moreover let S denote the solution set of (StQP):
S := arg min x
©1
2x
T
P x + c T x + Q(x) : Ax = b, x ≥ 0.
ª
We state the following result and refer the reader to [81, 82] for proof. The proof in [81, 82] uses the recourse function defined in terms of the solution the Dorn dual [26] of (subω ), given the strong duality between (subω ) and its Dorn dual. 30
Although we have used the solution of the Lagrange dual to define Q(x), the result from [81] continues to hold since strong duality applies in our case too. Proposition 1.4.2 The random recourse function Q(x; ω) is finite for all x ∈ K 2 (ω). Moreover Q(x; ω) is convex and lower semicontinous for all x ∈ K 2 (ω) and for all ω ∈ Ω. For a finite Ω, the recourse function Q(x) defined as
Q(x) := IEω [Q(x; ω)] , ∀ x ∈ K 2 is convex and lower semicontinuous in x.
1.4.1 L-shaped Method with Inexact Cuts The L-shaped method creates an outer-approximation of the recourse function by using a series of optimality cuts [13, 78]. The usage of inexact optimality cuts to model the recourse function for stochastic linear programs was discussed separately by Au et al. [4] and Zakeri et al. [98]. This technique showed an improvement in the performance of the traditional L-shaped method [92] since it required the solution of the dual problems only to feasibility, as opposed to optimality. The L-shaped method with inexact cuts was extended to solve stochastic QPs by Shanbhag et al. [81, 82]. In such a method, the function Q(x) is approximated by inexact optimality cuts of reducing inexactness. In the event that a choice of the first stage variable renders any (subω ) infeasible (or a corresponding dual problem (dsubω ) unbounded), the algorithm adds a feasibility cut. The assumption 1.3.1 of relatively complete recourse for (SNLP) does not imply that the QP approximation (StQPk ) would have relative complete recourse. This necessitates that infeasibility in a subproblem be handled by suitably restricting the first stage decisions through feasibility cuts. Next, we define an ²-inexact optimality cut and a feasibility cut based on [81] but using the solutions to the Lagrange duals (dsubω ). Definition 1.4.3 Suppose each dual problem (dsubω ) is solved to an optimality tolerance of ², with solution given by (v ω∗ , u ω∗ ). Then an ²-inexact cut is defined by G IT x + g I + ² > Q(x),
31
where G I and g I are given by £ ¤ IEω A Tω v ω∗ −1 ω∗ −1 T ω∗ −1 −1 T ω∗ u − 21 v ω∗T B D ω B v + v ω∗T B D ω d ω + u ω∗T D ω dω IEω [b ω v − 12 u ω∗T D ω −1 T ω∗ − u ω∗ D ω B v ],
respectively. Definition 1.4.4 Given an xe that solves (master), suppose the scenario dual (dsubω )
be unbounded. If (v¯ ω , u¯ ω ) is a direction of unboundedness for (dsubω ), boundedness of (dsubω ) is ensured by the following inequality −1 −1 −1 T ω d ω + (u¯ ω )T D ω d ω − u¯ ω D ω B v¯ ≤ 0. (b ω − A ω xe)T v¯ ω + (v¯ ω )T B D ω
A feasibility cut is thus defined as F ωT x ≥ f ω where T ω −1 −1 −1 T ω F ω := (A Tω v¯ ω )T and f ω := b ω v¯ + (v¯ ω )T B D ω d ω + (u¯ ω )T D ω d ω − u¯ ω D ω B v¯ .
Let I opt and I fea be the set of optimality and feasibility cuts present at the current iterate of the algorithm. Note that I fea contains elements identified by a tuple ( j , ω) consisting of an index and the scenario number. Note that if the directions of unboundedness form a multidimensional cone, then multiple feasibility cuts can arise out of the unboundedness of the same (dsubω ). In effect, the set I fea may have tuples with multiple indices for the same ω. At any iteration, the stochastic QP master problem, together with the cuts is then given by m-ils
minimize x,θ I
1 T T I 2x Px +c x +θ
Ax = b, subject to
jT G I x + θI F jT,ω x
j
≥ gI ,
j ∈ I opt
≥ f j ,ω ,
j , ω ∈ I fea
x ≥ 0. In the L-shaped method with inexact cuts (ILS), all the dual problems are solved to ²-inexactness to obtain the cut (G I , g I ). The inexactness is steadily reduced through successive iterations in such a way that {²k } that converges to zero [81]. This leads to the following definitions of upper and lower bounds. 32
Definition 1.4.5 The lower and upper bounds at iteration k are defined as 1 1 I L kI ≡ c T x k + x kT P x k + θk and UkI ≡ min{Uk−1 , c T x k + x kT P x k + Q(x k )}, respectively. 2 2 Notice that {L kI } is a monotonically increasing sequence while {UkI } is a monotonically decreasing sequence. The L-shaped method with inexact cuts (Algorithm 2) proceeds as given below [81]. Algorithm 2: L–shaped method with inexact cuts 0
initialization k = 1,UkI = ∞, L kI = −∞; choose ²1 , τ, u > 1; while |UkI − L kI | > τ do
1
Solve (m-ils) to get (x k , θk );
2
Update lower bound L kI ;
3 4
Solve (dsubω ) to tolerance ²k to obtain (u ω∗ , v ω∗ ) for all ω ∈ Ω; if (dsubω ) is bounded ∀ω then ¢ ¡ Construct G Ik , g Ik as per definition 1.4.3; else Add feasibility cut and go to step 1; end
5
¡ ¢ Update upper bound UkI and add optimality cut G Ik , g Ik to (m-ils);
6
Update inexactness ²k+1 = ²k /u;
7
k = k +1 end
In the above algorithm, u is called the inexactness update. It controls the rate at which optimality cuts are tightened. It can be shown that as k → ∞, x k approaches x∗ , the solution of (StQP) and θ approaches Q(x∗ ). The convergence k
theory for the inexact cut L-shaped method can be found in [81, 82].
1.4.2 Trust-region L-shaped Method for Stochastic QPs In [56], a trust-region variant of the L-shaped method is employed for solving stochastic linear programs. Specifically, such a scheme requires appending a trust-region bound to the master proble to regulate the step taken. In the past, trust-region approaches have been explored by Kiwiel [50] (bundle methods) and by Hiriart-Urruty and Lemaréchal [40]. 33
The trust-regions used by Linderoth and Wright differ from these in the way that they are box constraints, keeping the new set of constraints also linear. Furthermore, they assume that their stochastic linear program has complete recourse [56]. As noted in the previous section, assuming relatively complete recourse for (SNLP) does not guarantee us complete recourse for (StQPk ). Here we present a modified version of their algorithm for stochastic QPs that does not assume complete recourse. Our method has the following distinctions from the method of LW. • It solves stochastic QPs. • It does not assume relatively complete recourse and uses feasibility cuts Consider a generic stochastic quadratic program (StQP). If we assume that all (subω ) are feasible, solving (StQP) is equivalent to solving the following problem. master
minimize x,θ
1 T T 2x Px +c x +θ
Ax = b subject to
θ ≥ IE
£1 2
y¯ωT D ω y¯ω + d ωT y¯ω
¤
x ≥ 0. This is an exact-cut analogue of (m-ils) with y¯ω as the solution to (subω ). The right hand side of the second constraint forms the recourse function defined by definition 1.4.1 and theorem 1.4.2. One of our key modifications to LW’s method is that in the event that any (dsubω ) turns out to be unbounded, (i.e. (subω ) is infeasible) and y¯ω does not exist, we append the master problem with a feasibility cut as in [81], the form stated in definition 1.4.4. It is hence possible to append the feasibility cut to the constraint Ax = b in the form of an equality constraint with slack variables to form a new set of constrainsts ¯ In lemma 1.4.11, we show that the total number of feasibility cuts ¯ − s = b. Ax that can be added in the algorithm is finite. For simplicity, we assume that all the possible feasibility cuts are added and included within the constraint ¯ ¯ − s = b. Ax Before proceeding further, we introduce some notation. Let 1 b Q(x) ≡ x T P x + c T x + Q(x), 2 34
¡ ¢ where Q(x) is defined in definition 1.4.1. Let G k , g k be the k t h set of optimal-
ity cuts obtained from exact solution of all subproblems (dsubω ). As in [56], we define two kind of iterates – major and minor, of which the major fall into two types to be defined later. The sets I opt and I fea again represent the set of optimality and feasibility cuts respectively while S is the set of all solutions of (StQP). The model function formed by the set of optimality cuts that approxib mates Q(x) at the k t h major iteration is denoted as m k (x). Its domain is n .
R
m k (x) ≡ =
1 T x P x + cT x + 2
© ª max −G kT x + g k k∈I opt © ª T T opt 1 T x P x + c x + inf θ : θ ≥ −G x + g ∀ k ∈ I . k k 2 θ
(1.11)
Let ∆k denote the trust-region imposed at the k t h major iteration. Then, the first-stage problem at any major iteration, together with the box-shaped trust region is given by (m-trk ): m-trk
minimize x,θ
1 T x P x + cT x + θ 2
¯ − s = b¯ Ax subject to
θ ≥ −G Tj x + g j
∀ j ∈ I opt
−∆k e ≤ x − x k ≤ ∆k e x, s ≥ 0.
An alternative to using the `∞ trust regions would be to use `2 trust regions by adding a constraint to (m-trk ) like kx − x k k2 ≤ ∆k . While such a quadratically constrained QP would not be solvable by standard QP solvers, this is more in line with conventional trust-region subproblems [68]. Another alternative is to impose the trust region implicitly by adding a term ρkx − x k k2 as a penalty to the objective of (m-trk ). This results in a problem similar to the one in the regularized decomposition method [78]. We leave the question of which trust-region framework may prove advantageous as future research and focus primarily on extending the approach suggested in [56]. Let x¯ be the solution of (m-trk ) or © ª ¯ −∆k e ≤ x − x k−1 ≤ ∆k e, x, s ≥ 0 . ¯ − s = b, x¯ = arg min m k−1 (x) : Ax x
35
(1.12)
Then x¯ is taken to be the new major iterate of type 1 if and only if the acceptance test stated below is satisfied. We denote the k t h major iterate by x k , regardless of its type. Solutions x¯ that are not major iterates are called minor iterates. The l t h minor iterate following the k t h major iterate is denoted by x k,l , l = 0, 1, 2, . . .. The model functions formed in this sequence are denoted by m k,l (x). In the next subsection, we elaborate on how the model function gets updated particularly if the acceptance test is not met. Model update procedure We begin by noting that m k,0 (x) := m k (x). But note that a major iterate x k is not the 0th minor iterate x k,0 . In fact, x k is a solution to the model function given by m k−1,r where r represents the terminating iterate of the previous subsequence. The new model function adds the optimality cut and its solution is x k,0 . If for any first-stage solution x k,l , the dual problem (dsubω ) is unbounded for some ω ∈ Ω, we add a feasibility cut to (m-trk,l ) and form (m-trk,l +1 ). As a consequence, the sequence of minor iterates is taken to have terminated at x k,l −1 . The solution of (m-trk,l +1 ) is still denoted as x k,l since it minimizes m k,l (x), although the minimization is subject to the added feasibility cut. We continue to rename x k,l in this manner until for some value of x k,l , all the subproblems (dsubω ) are bounded for all ω ∈ Ω. This is taken as the major iterate of type 2. In effect, we terminate the subsequence when no further feasibility cuts are added. This is formalized in the following definition. Definition 1.4.6 Let x k,l be the current iterate and suppose for some ω ∈ Ω (dsubω ) 1 2 parametrized by x k,l is unbounded. Let x k,l , x k,l , . . . denote the solutions of subsequent problems formed by the addition of feasibility cuts.3 Suppose i¯ is the ¯
i i is taken to be a masmallest i such that all (dsubω (x k,l )) are feasible. Then x k,l
jor iterate of type 2 and renamed as x k+1 . Furthermore the sequence of minor iterates is said to have terminated at x k,l −1 . Note that a major iterate of type 2 need not satisfy the acceptance test. Lemma 1.4.11 shows that the number of feasibility cuts that can be added during the algorithm is finite. So a major iterate of type 2 is a well defined entity. To summarize, an iterate x k could be a major iterate of either type 1 or type 2. If the iterate has satisfied the acceptance test, then it is a type 1 iterate; if it is an 3
These values are stored in the variable x k,l .
36
iterate obtained by solving a master problem with a newly added feasibility cut for which all subproblems are feasible, then it is a type 2 iterate. We note that LW have assumed ([56] assumption 1) that their problem has complete recourse and hence do need require feasibility cuts. The minor iterate x k,bl solves a similar problem (m-trk,l ), l ≥ b l given by n o ¯ kx − x k k∞ ≤ ∆k,l , x, s ≥ 0 . ¯ − s = b, x k,bl = arg min m k,bl (x) s.t. Ax x
(1.13)
As per the above description, b l < l iff (m-trk,l ),(m-trk,l −1 ),. . ., (m-trk,bl ) differ only by feasibility cuts. The following procedure is similar to that in LW [56] with modifications made to accomodate feasibility cuts. To obtain m k,l +1 from m k,l , we proceed as follows: 1. If any scenario dual (dsubω ) is unbounded, then a feasibility cut is added to (m-trk,l ) and no cuts are deleted. The model function m k,l (x) remains unchanged as no optimality cuts are added. 2. If all (dsubω ) are bounded, (a) all cuts generated at x k are retained. (b) all cuts active at the solution of (m-trk,l ) are retained. (c) for a given η ∈ [ξ, 1), all cuts generated at previous minor iterates, lˆ during the same major iteration, that satisfy i h b k ) − m ˆ(x ˆ) m k,l (x k ) − m k,l (x k,l ) > η Q(x k,l k,l
(1.14)
are retained. (d) if a cut does not satisfy (a), (b) or (c) and has been inactive for the past 100 solutions of (m-trk,l ), then it is deleted. (e) all optimality cuts arising from the solutions of (dsubω ) at x k,l are added to m k,l (x) to get m k,l +1 (x). Note that since all cuts at x k (either a type 1 or type 2 major iterate) are retained, we have b k ), m k,l (x k ) = Q(x
37
∀ k, l .
(1.15)
To elaborate further, at x k the optimality cuts form an exact estimate of the recourse function. Thus the equality. Since the model keeps all the cuts at x k , this equality is maintained throughout the subsequence. Specifically b k ). m k,0 (x k ) = Q(x
(1.16)
Note from (1.13) that x k,l minimizes the model m k,l (x). If we consider the subsequence of iterates x k,l beginning with a major iterate x k , we have that the model function value at x k,l , namely m k,l (x k,l ), cannot be equal to the value of the recourse function value at that point or b k,l ) m k,l (x k,l ) 6= Q(x
x k,l ∉ S.
In fact, since a cut at x k,l is added to form the model function m k,l +1 (x), this b model is an exact estimate of Q(x) at x k,l . That is to say, b k,l ). m k,l +1 (x k,l ) = Q(x b b Since Q(x) is convex, m k,l (x) is a piecewise linear underestimate of Q(x). The
algorithm needs an initial model function m 0,0 to begin the algorithm. We require it to satisfy the following assumptions. b 0 ) where Assumption 1.4.7 At the initial starting point x 0 we have m 0 (x 0 ) = Q(x b m 0 (x) is a piecewise linear underestimate of Q(x)
Next, we define the acceptance test which dictates one avenue for the termination of the subsequence of iterates x k,l . Acceptance test For some fraction ξ ∈ (0, 1/2), the minor iterate x k,l is accepted if and only if b k ) − Q(x b k,l ) Q(x ¡
m k,l (x k ) − m k,l (x k,l )
¡ ¢ b k,l ) ≤ Q(x b k ) − ξ m k,l (x k ) − m k,l (x k,l ) . (1.17) ¢ ≥ ξ or Q(x
The trust-region approach limits the step size and thus allows us to solve (m-tr) with a limited number of optimality cuts. As a result, cuts that are inconsequential for the solution of x k,l can be deleted. Which cuts to delete is a matter of careful choice since incorrect deletion can make us re-evaluate cuts that were 38
deleted prematurely. For iterates that do not satisfy the acceptance test, LW provides a prodecure to update the trust region which we adopt and describe next. Updating the trust-region ∆ Let ρ be defined as ρ := min(1, ∆k,l )
b k,l ) − Q(x b k) Q(x . b k ) − m k,l (x k,l ) Q(x
(1.18)
Then the update of the trust region is carried out as follows: (1.) If [ρ > 0 for 3 iterations with same ∆k,l and ρ ∈ (1, 3]] or [ρ > 3], then the updated trust region is given by ∆k,l +1 =
1 ∆k,l . min(ρ, 4)
(2.) When the new major iterate x k+1 is obtained, there is a case for increasing the trust region. If ¡ ¢ b k,l ) ≤ Q(x b k ) − 0.5 Q(x b k ) − m k,l (x k,l ) and kx k − xk∞ = ∆k,l , Q(x ¡ ¢ then ∆k+1,0 = min ∆hi , 2∆k,l where ∆hi is a pre-specified upper bound
on the trust region. (3.) If (1.) or (2.) are not satisfied, then ∆k,l is left unchanged Finally, given ²t ol > 0 the algorithm terminates if ¡ ¢ b k ) − m k,l (x k,l ) ≤ ²t ol 1 + |Q(x b k )| . Q(x
(1.19)
We now state the complete L-shaped trust-region method (Algorithm 3) for solving stochastic QPs.
39
Algorithm 3: Trust-Region L–shaped method 0 initialization k = 1, l = 0; choose ξ ∈ (0, 1/2), initial point x 0,0 , initial model m 0,0 and initial trust region ∆0,0 ∈ (0, ∆hi ], ²t ol < 1; b k ) − m k,l (x k,l ) > ²t ol (1 + |Q(x b k )|) do while Q(x 1
Solve (1.13) to get x k,l ;
2
¡ ∗ ∗¢ Solve (dsubω ) accurately to obtain solution points u ω , v ω for all
ω ∈ Ω; 3
if (dsubω ) is bounded ∀ω then ¡ ¢ Construct G k , g k ; if feasibility cut was added in the previous iteration then x k+1 = x k,l ; Obtain new model m k+1,0 as per section 1.4.2; end else Add feasibility cut and go to step 1; end
4
if acceptance test (1.17) is satisfied then Set x k+1 = x k,l ; Obtain the new model m k+1,0 as per section 1.4.2; Update ∆k+1,0 as per section 1.4.2; else Update the model to get m k,l +1 as per section 1.4.2; Update ∆k,l +1 as per section 1.4.2;
5
end 6
k = k + 1; end
Convergence of the trust region method An outline of the main convergence results follows. Since, this theory follows portions of that in LW closely, we carefully distinguish our results from theirs particularly when the proofs either differ or have been extended significantly. Our first two results establish that the model function does indeed increase with every minor iterate and provide a bound for such an increase. Next, we prove a boundedness result for the trust-region radius as well as the finite termination of the subsequence of iterates originating from a major iterate. Re-
40
lying on properties associated with the aforementioned subsequence, the convergence theorems for the algorithm are stated and proved. The model function is an under-approximation of cost of first and second-stage decisions (recourse). Since the model function is convex, the algorithm converges when the difference between the value of the model function at its minimum (m k,l (x k,l )) b k )) is within a pre-specified and its value at the last major iterate (m k,l (x k ) = Q(x tolerance level. Lemma 1.4.8 shows that the optimal value of the model function increases with every minor iterate. Lemma 1.4.8 If a minor iterate x k,l +1 does not satisfy the acceptance test (1.17), we have m k,l (x k,l ) ≤ m k,l +1 (x k,l +1 ). Proof : See lemma 1 in LW [56]. In lemma 1.4.9, we provide a lower bound on the difference between the optimal value of the model function and the value of the model function at the last major iterate, namely m k,l (x k ) − m k,l (x k,l ). Lemma 1.4.9 Consider an iterate x k,l derived at a major iteration k and minor iterate l where k, l ≥ 0. Then the following bound on the change in the model function holds: ∆k,l
µ b k ) − m k,l (x k,l ) ≥ min m k,l (x k ) − m k,l (x k,l ) = Q(x
kx k − P (x k )k∞
¶ ¡ ¢ b k ) − Qb∗ , 1 Q(x
where P (x k ) is the projection of x k on S. Proof : We begin by noting that x k,l minimizes m k,l (x) in problem (m-trk,l ). Hence for all x that are feasible for (m-trk,l ), we have m k,l (x k,l ) ≤ m k,l (x).
(1.20)
In particular, this holds for x = x k + ∆k,l
P (x k ) − x k . kP (x k ) − x k k
Recall that (m-trk,0 ) and (m-trk,l ) differ only in the optimality cuts. The addition of a feasibility cut creates a new major iterate. Hence x k is feasible for (m-trk,l ).
41
P (x k ) is a solution to the problem so it is feasible as well. In other words ¶ P (x k ) − x k ¯ A¯ x k + ∆k,l = b. kP (x k ) − x k k µ
From (1.20), it follows that ¡ ¢ m k,l (x k,l ) ≤ m k,l x k + βk (P (x k ) − x k ) where βk =
∆k,l kP (x k ) − x k k
.
b Since m k,l (x) is a lower underestimate of Q(x), it follows from the convexity of Qb that ¡ ¢ m k,l (x k,l ) ≤ Qb (1 − βk )x k + βk P (x k ) , b k ) + βk Q(P b (x k )) = Q(x b k ) + βk (Qb∗ − Q(x b k )), ≤ (1 − βk )Q(x
where b (x k )) = Qb∗ . Q(P
Using (1.16) and rearranging the above inequality, we obtain b k ) − m k,l (x k,l ) ≥ Q(x
∆k,l kP (x k ) − x k k
¡ ¢ b k ) − Qb∗ . Q(x
(1.21)
b k ) − m k,l (x k,l ) ≥ Q(x b k ) − Qb∗ . But kP (x k ) − x k k ≤ ∆k,l , implying that Q(x
Combining this with (1.21), we obtain the required result. Next, we prove that the trust-region is bounded away from 0 for all iterates x k,l ∉ S. Lemma 1.4.10 Suppose E k , F k and β are defined as E k ≡ min kx k¯ − P (x k¯ )k∞ ,
Fk ≡
¯ k=0,1,...,k
min
¯ k=0,1,...,k,x k¯ ∉S
b ¯ ) − Qb∗ Q(x k
kx k¯ − P (x k¯ )k∞
,
b β = sup kg k, g ∈ ∂Q(x), x feasible for (m-trk,l )
˜ l˜ such that min(E ˜ , F ˜/β) < ∆hi . Then for all k > k˜ and for Suppose there exists k, k l all l , the trust-region ∆k,l satisfies ∆k,l ≥ (1/4) min(E k , F k /β). Proof : See lemma 3 in LW [56]. 42
Lemma 1.4.11 is a result distinct from that proved by LW. It shows that the algorithm can have only finitely many feasibility cuts. Lemma 1.4.11 The subsequence of major iterates of type 2 terminates finitely. Proof : Recall that the addition of a feasibility cut relies on one of the subproblems, namely (subω ) being infeasible. Equivalently, this amounts to the corresponding dual problem (dsubω ) being unbounded. The region of unboundedness of each (dsubω ) is spanned by the set of extreme rays. Since the dimensionality of (dsubω ) is finite, the number of extreme rays is clearly finite as well. Each feasibility cut corresponds to one and only one extreme ray. Thus the maximum number of feasibility cuts is therefore finite by virtue of the finiteness of the cardinality of the sample space (|Ω| < ∞). It follows that the subsequence of major iterates of type 2 cannot be infinite. The next set of results require specifying the tolerance ²t ol in the algorithm for which we consider two cases: (1) ²t ol = 0: Here, we show that the sequence of major iterates either terminates finitely in a point in S or converges to a point in S. (2) ²t ol > 0: In this instance, we show that the algorithm terminates finitely to a point in S. Lemma 1.4.12 proves that the sequence of minor iterates either produces an iterate that satisfies the acceptance test or produces a minor iterate that results in a strict reduction in m k,l (x k ) − m k,l (x k,l ). Lemma 1.4.12 Let ²t ol ≥ 0 and ξ be given where ξ specifies the relative improvement in the acceptance test in (1.17). Let η be the parameter used for specifying the deletion of cuts as per (1.14). For some major iterate k of type 1 and minor iterate l 1 , suppose x k,l 1 fails the acceptance test given by (1.17). Then the sequence of minor iterates generates an x k,l 2 that satisfies (1.17) or there is an index l 2 that satisfies £ ¤ b k ) − m k,l (x k,l ) ≤ η Q(x b k ) − m k,l (x k,l ) . Q(x 2 2 1 1
Proof : We consider two cases: (a) ²t ol = 0: See lemma 4 in LW [56] .
43
(1.22)
(b) ²t ol > 0: This result holds by observing that the relation ζ :=
η−ξ b k ) − m k,l (x k,l )] > 0, [Q(x β
in expression (52) of LW [56] holds even when ²t ol > 0. For the final convergence statement, we first consider the case ²t ol = 0. Combining the implications of lemma 1.4.12 with lemma 1.4.9 we can show that, if x k ∉ S, then it produces a minor iterate x k,l that satisfies the acceptance test. This is the result of theorem 1.4.13. Theorem 1.4.13 Suppose ²t ol = 0. Then the following hold: (1) If x k ∉ S then there exists an l ≥ 0 such that x k,l satisfies (1.17) (2) If x k ∈ S then either the algorithm terminates finitely or ¢ ¡ b k ) − m k,l (x k,l ) Q(x
decreases monotonically to zero. Proof : See theorem 1 in LW [56]. Theorem 1.4.14 shows that the sequence of major iterates approaches a solution. Theorem 1.4.14 Suppose ²t ol = 0. The sequence of major iterations is either finite terminating at some x k ∈ S or infinite with kx k − P (x k )k∞ → 0. Proof : See theorem 2 in LW [56]. We conclude our convergence results with a finite termination result for the case where ²t ol > 0. This result was claimed in LW to be have been easy to prove and skipped. But since our algorithm differs from LW’s to a certain extent we prove this result here for completeness. Theorem 1.4.15 For ²t ol > 0, the algorithm terminates finitely with a solution to (StQP). Proof : Assume in accordance with lemma 1.4.11 that all feasibility cuts have been added. From lemma 1.4.12, for some minor iterate l that fails the acceptance test (1.17), the sequence of minor iterates either yields an iterate that satisfies test (1.17) or an iterate satisfying Eq (1.22). Let x k,l 1 be the current iterate that does not satisfy test (1.17). Two cases arise: 44
(a) Suppose that the sequence of minor iterates emerging from major iterate k never satisfies the acceptance test. Then for l 2 > l 3 > . . . l r we have ¡ ¢ ¡ ¢ b k ) − m k,l (x k,l ) ≤ η Q(x b k ) − m k,l (x k,l ) ≤ . . . ≤ ηr −1 Q(x b k ) − m k,l (x k,l ) . Q(x r r r −1 r −1 1 1 b k) − This indicates that for sufficiently large, but finite r , the term Q(x
m k,l r (x k,l r ) can be reduced enough to satisfy (1.19), implying the finite termination of the algorithm. (b) Suppose the test (1.17) is satisfied for some l 1 . This results in a new major iterate x k+1 . If for any k¯ ≥ k the sequence of minor iterates x ¯ , l ≥ 0 k,l
behave as in case (a), then convergence follows from the argument in (a). So suppose that for all k¯ ≥ k the minor iterate sequence x ¯ , l ≥ 0 does not k,l
behave as in case (a). This means that for all k¯ ≥ k the acceptance test is satisfied for some x k,l ¯ , l ≥ 0. For major iterates that have not terminated we have b ¯ ) − m ¯ (x ¯ ) > ²t ol (1 + |Q(x b ¯ )|) Q(x k k,l k,l k
Hence every major iterate causes a reduction in Qb of at least ξ²t ol as seen below ´ ³ b ¯ ) − m ¯ (x ¯ ) ≥ ξ²t ol . b ¯ ) − Q(x b ¯ ) ≥ ξ Q(x Q(x k k,l k,l k k,l b ¯ ) is bounded Since a solution exists to the minimization problem (StQP), Q(x k
below. This further implies that the minimum of (StQP) is achieved in finitely many major iterations.
1.4.3 Recovering the Multipliers A drawback of using the L-shaped method to solve (StQP) is that the multipliers for (StQP) may not be directly obtained. The master problem and the subproblems are solved by QP solvers which return multipliers to their respective constraints, which are not the multipliers to the original problem. To explain this, we begin with some notation. All variables with an asterix as superscript pertain to the KKT point of (StQP). Hence we have the conditions for optimality
45
of (StQP) as P x ∗ + c − A T λx∗ −
X ω
A Tω λω∗ − µx∗ = 0
£ ¤ p ω D ω y ω∗ + d ω − B T λω∗ − µω∗ = 0 ∀ ω ∈ Ω
Ax ∗ = b B y ω∗ = b ω − A ω x ∗
∀ω∈Ω
0 ≤ µx∗ ⊥ x ∗ ≥ 0 0 ≤ µω∗ ⊥ y ω∗ ≥ 0 ∀ ω ∈ Ω. The solutions of the subproblems (dsubω ) are denoted by a superscript s, and at a given x they satisfy their respective KKT conditions £ ¤ p ω D ω y ωs + d ω − B T λωs − µωs = 0
B y ωs = b ω − A ω x 0 ≤ µωs ⊥ y ωs ≥ 0,
∀ ω ∈ Ω.
The master problem solution (for the inexact L-shaped method) is denoted by a superscript m. Assume for simplicity that there are only optimality cuts. If 1 represents the column of ones, then the KKT point of the master satisfies P xm + c −
X j
T xm G j γm =0 j −A λ
1−
X j
γm j =0
Ax m = b jT m 0 ≤ γm x +θ−g j ≥ 0 j ⊥G
0 ≤ µm ⊥ x ≥ 0. At the optimal point we have x m = x ∗ , but the multiplier λxm satisfies a different set of equations from λx∗ and hence λxm 6= λx∗ and µxm 6= µx∗ . However, since the linear independence constraint qualification holds for (StQP) (in this case requiring that B ω has full row rank for each ω ∈ Ω), when y ωs = y ω∗ , it follows that λωs = λω∗ and µωs = µω∗ , ∀ ω ∈ Ω. The outer loop of the SQP method requires that we have the accurate value of λx∗ . This section is devoted to recovering the accurate value of the λx∗ . Let
46
z,Q, q, M , h, λ and µ be defined as
x
z :=
y1 .. .
M
y |Ω|
P 0 Q := 0 . ..
A
0
:=
A1 .. .
B
A |Ω|
0
0
...
0
T c p1D 1 ... 0 p 1 d1 0 p2D 2 ...0 .. q := . .. .. . . p |Ω| d |Ω| 0 ... 0 p |Ω| D |Ω| ... 0 b λx µx 1 1 b1 λ µ ... 0 h := λ := µ := . , respectively. .. .. .. . . 0 . . . ... B b |Ω| λ|Ω| µ|Ω|
It follows that (StQP) is given by StQP
minimize z
1 T T 2 z Qz + q z
subject to
Let
Ã
N=
MQ −1 M T
Q −1 M T
Q −1 M T
Q −1
Mz = h
:λ
z ≥0
:µ
!
Ã
η=
MQ −1 q + h
!
and ν =
Q −1 q
à ! λ
µ
.
Then the Lagrange dual of (StQP) is StQPdual
maximize − 12 νT N ν + ηT ν − 12 q T Q −1 q λ,µ
subject to
µ ≥ 0,
P where p = dim(b)+ ω dim(b ω ). At z = z ∗ , solution to (StQPdual) is λ = (λx∗ , λ1∗ , . . . , λ|Ω|∗ )
and µ∗ = (µx∗ , µ1∗ , . . . , µ|Ω|∗ ). Now from the solution of the L-shaped method, we know the accurate values of λω∗ , x ∗ , y ω∗ . Thus this problem boils down to a quadratic in only λx , µx . Ã MULT
maximize x x λ ,µ
subject to
λx
!T
µx
e N
µx ≥ 0
47
à ! λx
µx
à e +η
T
λx µx
!
where à e= N
AP −1 A T
P −1 A T
P −1 A T
P −1
!
à e= and η
−AP −1
P
A Tω λω + b + AP −1 c − P −1 P −1 c
P
A ω µω
!
.
It is noteworthy that this involves solving only one more QP for every SQP iteration and hence does not significantly impact the overall performance of the algorithm.
1.5 Numerical Results This section discusses the numerical performance of our algorithm from a variety of standpoints. We begin by examining the behavior of the algorithm with both subproblem solvers on a set of stochastic nonlinear programming test problems. Next, we discuss the impact of changing key parameters associated with the workings of the subproblem solvers. Finally, this is followed by a comparative study between the performances of the inexact cut and the trustregion variants of our proposed algorithm.
1.5.1 Scalability of the Algorithm There is no standard stochastic NLP test problem set available. Chen and Womersley [18] have developed test problems for stochastic quadratic programs. We use this set to generate two sets of test problems (convex and general nonlinear) as follows. • Convex problems: The objective function was taken to be the norm function (of 2nd and higher order), parametrized by the scenario number ω. The constraints were taken from the Chen and Womersley set for the corresponding number of realizations. This way we ensure that we have the required recourse properties on the subproblems. Specifically, the objective function is given by f (x) = 12 x T H x + c T x + d + (1 +
X (x(i ) − a i (ω)))1/p
where H Â 0, p = 2, 4, . . . and a i (ω) is a function of ω.
48
• Nonlinear problems: In this case we used the first stage objective as a function from the Hock–Schittkowski set [44] and the second stage objective as a convex function. The constraints were again taken from the Chen and Womersly set. Following are the two problems we used along with their index in the Hock–Schittkowski set. f (x) = (x(1) − x(2))2 + (x(3) − 1)2 + (x(4) − 1)4 + (x(5) − 1)6
(Prob.35)
f (x) = (x(1) − 1)2 + (x(1) − x(2))2 + (x(2) − x(3))2 + (x(3) − x(4))4 + (x(4) − x(5))4
(Prob.58).
In the results that follow, we present some representative problems focusing mainly on the effect of changing |Ω| and the comparative behavior of the two algorithms. cpu denotes the CPU time in seconds and maj denotes the number of outer SQP iterations. The other fields are self-explanatory. For every case, the first stage variable x and the second stage variable y ω had dimension 5. The constraints were linear, with 4 rows for the first stage constraint and 5 rows for the second stage constraints. Thus the number of variables = 5 × |Ω| + 5. The number of constraints = 5 × |Ω| + 4, in addition to 5 × |Ω| + 5 nonnegativity constraints. Unless otherwise stated in the tables, all algorithmic parameters (like inexactness update, ∆hi etc.) were held constant across problems. The convex problems were solved using exact hessians, while the nonlinear problems used our update rule. In most cases, the KKT residual desired to satisfy (1.10) is about 10−5 . The reader may observe that for several instances the KKT residual obtained has been greater than this number. This is because the algorithm has been terminated because of minimal progress. A detailed assessment of this phenomenon will follow later. It is notable that both TR and ILS algorithms show good scaling with |Ω|. We have noted some problems to show have taken far less effort than expected which we can regard as being exceptions from the broad linearity seen otherwise. From our tables it can be seen that for convex problem ½
max
qp1 /qp2
¾
|Ω1 |/|Ω2 |
≈
2,
ILS
3,
TR.
Expectedly the performance on nonlinear problems is worse, though still linear. We also observed that for convex problems, the both solvers showed little 49
Table 1.1: Computational effort for the ILS solver – with convex costs 2, a i (ω) = ω/|Ω| ∀ i , u = 2) Prob 1 2 3 4 5 6 7 8 9 10
|Ω| 729 1000 1225 1600 2187 6561 16384 32768 59049 65536
cpuils
majils
2432 1446 3793 1293 3643 14967 74734 63300 218836 64599
6 5 4 3 5 5 5 4 4 3
cutsils
qpsils
kktils
348 150 321 83 174 217 430 189 354 87
254046 150155 393229 132886 380717 1423959 7045555 6193345 20903704 5701722
2.2E-04 6.0E-05 1.0E-06 2.7E-05 1.9E-04 3.8E-04 9.6E-05 1.1E-04 2.1E-05 3.0E-06
Table 1.2: Computational effort for TR solver – with convex costs 2, a i (ω) = ω/|Ω| ∀ i , ∆hi = 10) Prob 1 2 3 4 5 6 7 8 9 10
|Ω| 729 1000 1225 1600 2187 6561 16384 32768 59049 65536
cputr
majtr
2676 1697 13868 1026 4368 11141 35145 27764 285652 51568
6 8 12 3 4 6 3 3 6 3
cutstr
qpstr
kkttr
377 172 1165 65 208 174 223 81 482 70
275216 172180 1428302 104068 455108 1141794 3653858 2654292 28462106 4587593
1.0E-08 1.5E-04 2.4E-05 2.4E-04 1.6E-04 1.0E-06 7.4E-05 1.1E-04 4.8E-05 1.3E-04
(p =
(p =
change in the number of major SQP iterations. For nonlinear problems, we found no such trend. An important insight obtained from the results was that TR solver provided better KKT points (lower residuals) than the ILS method, but the ILS method converged faster.
1.5.2 Modifying the TR Update and Inexactness Update An interesting trend emerges for the trust-region method when the maximum trust region ∆hi is varied. Recall that at any iteration we must have ∆k,l ≤ ∆hi . Let us denote ∆∞ = max{∆k,l } k,l
as the largest trust region used by the algorithm during an implementation with ∆hi = ∞. We carried out tests for a fixed convex problem by increasing ∆hi from 5 to 500. The results are tabulated in table 1.5.2. It was seen that an increase in ∆hi lead to greater number of cuts and longer 50
Table 1.3: Computational effort for the ILS solver – nonlinear costs u = 2) Prob 1 2 3 4 5
|Ω| 729 1000 1225 1600 2187
cpuils 7336 3842 5266 5435 10844
majils 23 24 13 29 26
cutsils
qpsils
kktils
1010 385 430 337 501
737323 385409 527193 539566 1096214
9.4E-03 2.4E-04 4.6E-05 9.2E-04 5.2E-04
Table 1.4: Computational effort for the TR solver – nonlinear costs ∆hi = 10 ) Prob 1 2 3 4 5
|Ω| 729 1000 1225 1600 2187
cputr
majtr
11254 4706 5549 12607 16843
31 27 23 50 39
cutstr
qpstr
kkttr
1589 474 460 791 787
1159970 474474 563960 1266391 1721956
5.0E-06 3.0E-03 1.5E-05 1.3E-04 1.3E-04
(Prob. 35,
(Prob. 35,
CPU times. This behaviour is expected. In the implementation of the trust region method, we begin with a trust region of size 10−3 . Such a small trust region often leaves the initial point infeasible. In such a situation we increase the trust region gradually until a feasible point is obtained. This trust region is ∆0,0 . After the first iterate the trust region is increased only when a new major iterate is found and the model update rule for this case is satisfied. If ∆hi is large, the model update rule permits a large ∆k,l leaving the algorithm with a large region to be approximated by cuts, leading to more cuts. After a certain value an increase in ∆hi does not make any difference to the performance since we get max{∆k,l } = ∆∞ ¿ ∆hi . k,l
Table 1.5: Variation of computational effort and accuracy with ∆hi for TR method Prob 1 2 3 4 5 6 7
|Ω| 729 729 729 729 730 729 729
∆hi 5 10 15 20 30 50 500
cpu
maj
cuts
qps
3157 2676 5947 9056 9013 8971 9097
6 6 8 10 10 10 10
240 377 456 685 685 685 685
175200 275216 332880 500050 500735 500050 500050
maxk,l {∆k,l } 5 10 15 16.38 16.38 16.38 16.38
kkt 1.0E-08 1.0E-08 2.0E-06 1.0E-06 1.0E-06 1.0E-06 1.0E-06
For the ILS method, we varied the inexactness update (u in the ILS algorithm, denoted by upd in the table below) for a convex problem with 1600 sce51
narios. The inexactness was varied from 1.2 to 100. The results are tabulated in table 1.5.2. While changes in the update from 1.2 to 4 made a difference to the performance of the algorithm, any further increase resulted in no change. Specifically, we see that the algorithmic performance worsened by a more aggressive tightening of cuts. More cuts were needed, leading to more QPs being solved and we also ended up with worse KKT residuals. In addition, for larger u, the cuts were closer to exact implying that one needed to dual optimal solutions, rather than just dual feasible solutions. Further increases in u makes little incremental difference to the number of near–exact cuts that the algorithm deals with. Table 1.6: Variation of computational effort and accuracy with inexactness update for ILS method Prob 1 2 3 4 5 6 7 8 9 10
|Ω| 1600 1600 1600 1600 1600 1600 1600 1600 1600 1600
upd
cpu
sqp
cuts
qps
kkt
1.2 1.5 2 2.5 3.5 4 8 10 15 100
4569 2229 1293 1744 1595 3664 3101 3215 3911 3179
3 3 3 5 5 5 5 5 5 5
293 137 83 110 100 100 100 100 100 100
469093 219337 132883 176110 160100 160100 160100 160100 160100 160100
1.1E-05 9.0E-06 2.7E-05 1.3E-04 3.0E-04 2.8E-04 2.7E-04 2.7E-04 2.7E-04 2.7E-04
(a) QPs solved
(b) KKT residual
Figure 1.1: Performance profiles comparing ILS and TR solvers
52
1.5.3 Comparing ILS and TR In this section, we present two types of performance profiles [25]: 1. with the performance metric as the number of QPs solved 2. with the performance metric as the KKT residual (note that the x−axis is plotted in logarithmic scale). Our termination criteria for the algorithm were very strict, with τP = τD = 10−6 in criteria (1.10), while the termination criteria for the inner stochastic QPs were identical for both methods: τ = 10−10 for the ILS method and ²t ol = 10−10 for TR. It is clear that while the two methods are comparable purely on the basis of the number of QPs solved, the TR algorithm finds KKT points with lower residuals. Some explanations can be offered for this observed behaviour. The “acceptance test” in the TR method (1.17) is essentially a sufficient descent criterion. Theorem 1.4.13 shows that the convergence of b k ) − m k,l (x k,l ) Q(x
is monotone. The ILS method on the other hand offers no monotonicity guarantee. This is the crucial difference between the methods. While iterates of the ILS method might hover around the optimal point, the TR method proceeds steadily towards it, often leading to a more accurate KKT point of (SNLP). In the practical implementation of the algorithm we found that the Benders-SQP method may terminate due to the step length α obtained from the linesearch becoming arbitrarily small. As a consequence, the algorithm cannot make any further progress. Our understanding is that the latter becomes a possibility when the stochastic QPs are not solved to desired accuracy. As a result the multipliers obtained by solving problem (MULT) are not compatible with the primal solution, and what results is not a KKT point of the stochastic QP. This error propagates and eventually the algorithm ends up at a point from where it can neither make any progress, nor satisfy criterion (1.10). In such a situation we declare the algorithm to have terminated. Indeed, we have observed that tightening the termination criteria for the stochastic QPs results in better (lower) KKT residuals for (SNLP).
53
1.6 Summary, Concluding Remarks and Future Directions This paper has presented a novel technique of solving an important class of problems – stochastic nonlinear programs. The suggested method satisfies the key requirements articulated in the introduction. In particular, the algorithm displays both global convergence and superlinear local convergence properties. It is observed that the resulting QP approximations are two-period stochastic programs with random recourse. We provide a characterization for the firststage feasible decisions that in turn result in bounded feasible second-stage solutions with probability one. Finally, under the more practical assumption of a finite sample-space, we show that no extra assumptions on the recourse matrices are required. The third requirement on our framework is that of scalability. This is met through a departure from active-set QP solvers for solving the resulting QP approximations. Instead, we suggest two variants of cutting-plane methods. The first of these is an inexact variant of the L-shaped method first suggested in the second author’s dissertation [81]. The second solver is an extension of an LP variant in [56]. We extend the latter to the QP regime under weaker assumptions of recourse and provide convergence statements under these assumptions. Both these methods allow us linear scaling with the size of probability support, as our numerical results show. They however come with a disadvantage of not directly providing us with the multipliers for the constraints. To resolve this problem, we suggest the solution of single low-dimensional QP to obtain such multipliers. The structure of such a QP follows from Lagrangian duality and does not hamper the overall performance of the algorithm. Numerical results have shown the algorithm to have good scalability and convergence properties and particularly the ability to solve problems in order of a few hundred thousand variables and constraints. We consider this work the first step in the development of algorithms for solving truly large-scale multiperiod nonlinear stochastic optimization problems. Many questions however still persist. First, while the subproblem solvers satisfy scalability requirements, fast local convergence cannot be guaranteed. There appears to be some hope in such a direction based on work in bundle methods [33]. Second, do immediate extensions of the proposed framework
54
to the multiperiod domain exist?. Our view is that Benders framework has provided much by way of solving multistage LPs [13, 12]. Through QP extensions to such approaches, we believe that large-scale multistage stochastic NLPs could be solvable by Benders-SQP techniques.
1.7 Acknowledgments We are grateful for comments made by Stephen Wright and Suvrajeet Sen during the Informs Annual Meeting, 2007. We also appreciate the comments and suggestions provided by Sven Leyffer during his visit to UIUC in fall 2007.
55
2
Stochastic Nash-Cournot Equilibrium Problems: Theoretical Properties ———————————–
2.1 Introduction Planning1 models have evolved significantly since the last century when the focus lay on single-player deterministic models. Such problems were captured to some extent by linear programs [21] and more generally by nonlinear programs [31, 68]. However, such models have proved insufficient in the face of increasing complexity, specifically originating from competition, stochasticity and dynamics. Since the seminal work by Dantzig [21] and Beale [8] on planning under uncertainty, the field of stochastic programming has made significant inroads into areas of multi-period linear programming and discrete optimization (see [13]). More recently, planning models have undergone an evolution in that decision makers need to be cognizant of the behavior of one’s competitors in making strategic decisions. Specifically, the question has shifted from mere optima to that of equilibria. In particular, our interest is in obtaining Nash equilibria; namely a set of agent decisions at which no agent would like to deviate given the other agent’s decisions. Spurred by Nash’s ideas in the 50s [66], Nash equilibrium models have proved effective in capturing interactions in a variety of noncooperative settings ranging from electricity markets [27], telecommunication networks [1, 2] and traffic equilibria [65]. A key shortcoming in the models employed for capturing competitive settings lies in the suppression of the inherent uncertainty in the parameters. This is far from surprising since the resulting equilibrium problems can be large and therefore challenging to solve, particularly because the problem grows linearly with the number of players as well as with the cardinality of the sample-space. While research in stochastic programming has led to methods for solving prob1
This work comprising of chapters 2 and 3 will soon be submitted to Operations Research.
56
lems in linear, integer and more recently in the nonlinear settings, the Nash equilibrium problem with agents solving stochastic optimization problems has been left largely ignored. In this paper, we address precisely such a gap by focusing on the two-period stochastic Nash-Cournot equilibrium problem. A variety of Nash-Cournot settings exist in practice. Our focus is on developing tools for settings in which one has to participate in two or multiple markets, each of which is coupled to the previous. One obvious formulation for capturing such dynamics is through the framework for EPECs or equilibrium programs with equilibrium constraints [55]. Such problems model the Nash equilibrium problem of competing in the first-period market (say the forward market) and requires that each agent is constrained by equilibrium in the second-period market (such as the real-time market). Many challenges immediately follow from such a formulation. First, guarantees of existence of such equilibria are hard to come by in general cases. Second, even if one were to provide settings where such equilibria did indeed exist, the availability of globally convergent solvers has proved to be elusive for the present. Finally, manifold uncertainty in system specification make such approaches even less attractive. This paper pursues a different direction. Instead of solving the original dynamic game and determining a sub-game perfect Nash equilibrium, it solves a static game, albeit a much larger one. Specifically, we consider a multi-player decision-making problem in which every player competes in the first-period market with the assumption that he can take recourse in the second-period by modifying decisions with a corresponding cost. We define such a solution concept, an S-adapted open loop equilibrium, a generalization of S-adapted open loop control. This concept is explained in section 2.2.1. Noteworthy is the fact that such a formulation does not entail any ‘sequential decision making’ on the part of the players. The second-period decisions of players are dependent on the first-period decisions, but are not taken with any extra information available. As a result, the second-period decision of a player is not a singleton, but a collection of decisions he would take under different scenarios in the future. This important distinction will made clearer in a formal manner in a section 2.2.1. The contribution of this paper lies in formalizing our model in rigorous game–theoretic terms, studying its properties and developing a scalable algorithm for computing equilibria for our model. Broadly speaking, we consider a setting where agents are faced with determining participation levels in a mar57
ket whose price function is denoted by a non-increasing price function p(X ) where X=
n X
xi .
i =1
In keeping with the tenets of stochastic programming, our intent is to determine a robust set of bidding decisions which minimizes the cost of first and second-stage recourse actions. One formulation of the resulting class of problems of interest is given below. Definition 2.1.1 Consider a Nash equilibrium problem with N agents. Let agent i solve the parameterized stochastic optimization problem below. Agenti (x −i , y −i )
maximize p(X )x i − f i (x i ) + IE[p ω (Y ω )y iω − h iω (y iω )] x i ,y i
subject to
xi ∈ K i y iω ∈ K iω (x i ).
∀ω ∈ Ω,
where we use y i := (y i1 , y i2 , . . . , y i|Ω| ), the functions f i and h iω are assumed to be strictly convex twice continuously differentiable functions and K i are assumed to be a closed and convex set and K iω (x) is a closed convex set-valued map. Then a stochastic Nash equilibrium is said to be given by {x i∗ , y i∗ }iN=1 where ∗ ∗ (x i∗ , y i∗ ) ∈ SOL(Agenti (x −i , y −i )),
i = 1, . . . , N ,
and SOL represents the solution set of the optimization problem. Using this formulation as a basis, we now present the principal goals of this chapter: 1. First, we cast the above model in a formal game–theoretic framework. 2. Second, we prove existence and uniqueness of the stochastic Nash equilibrium for the above game. 3. Third, under fairly general assumptions, we develop an algorithm for obtaining such equilibria. our framework satisfies the key requirements of global convergence, fast local convergence and scalability with respect to the size of the underlying probability measure. 4. Fourth, we demonstrate the workings of this algorithm on a set of stochastic complementarity test-problems. 58
Items 3. and 4. above are included in the following chapter. The remainder of this chapter is organized as follows. In section 2, we provide some motivation for such a study through two instances of stochastic Nash-Cournot equilibrium problems. In addition, an outline of the literature on the Nash-Cournot framework is provided along with a review of solution concepts for stochastic games. Section 3 provides a formulation of the deterministic Nash-Cournot equilibrium problem. Section 4 is devoted to the stochastic Nash-Cournot problem. A sequential Newton-type matrix splitting algorithm is presented in the next chapter, together with the numerical results.
2.2 Motivation and Background Nash equilibrium problems arise in a variety of settings with a whole host of complexities. Section 2.2.1 comprises of three subsections, each of which provides a motivating example for studying such a class of problems. The first pertains to the modeling of competitive decision-making in an electricity market. The second problem represents a game-theoretic extension to the canonical newsvendor problem while last example is that of trading in emissions rights.
2.2.1 Source Problems In this section we outline some examples to motivate the stochastic Nash-Cournot equilibrium problem. The intention of this section to depict how stochastic Nash-Cournot games arise in practical settings. The games are hence described in their respective normal forms only. A more thorough game–theoretic treatment of these games and the type of equilibria we seek is done in the following section. Two-period electricity markets Here we consider a two-settlement Pool-type wholesale electricity market. In such a system, power generators bid into a day-ahead market (DA) to meet the forecasted load. An independent system operator (ISO) dispatches bids in a most economical way to clear the market, subject to physical transmission constraints. The selected bids are paid the market clearing price in the day-ahead market. When actual load is realized in the next day, we allow generators to
59
adjust their bids to clear the real-time market (RT). The differences of generators’ bids in the day-ahead market and in the real-time market are settled by the real-time market clearing price. Such a two-settlement system is used in several markets in the US, including ISO-NE, NYISO and PJM, while CAISO and ERCOT are moving toward this direction [49]. Here we consider that power generators participate in a two-settlement system, and bid into the day-ahead market to deliver energy in a particular hour on the following day. We assume that bidders have complete information about the game to be played; that is, bidders know the number of players in the game, the strategy space of each player, and the profit function of each player; but have incomplete information about events in the future. In addition, we assume that bidders share a common knowledge on the probability space (Ω, F, P ) of the uncertainties in the real-time market. In such a market, generators bid all their generations into the Pool. The ISO collects the bids and dispatches electricity to meet the demand in the most economical way. By doing so, the ISO can be viewed as auctioning off transmission services, while also performing spatial arbitrage to eliminate any price differentials in the network not caused by transmission congestion rents. For generators, we model a Cournot-type game; namely, each generator assumes that its bids on generation quantity would not affect other players’ actions on their quantity bids. For simplicity let us assume a single node network. Assume that market prices are determined by a linear inverse demand function at each demand node. Suppose N is a finite set of players. Let s i be the supply bid of firm i in the DA market. Let u iω , v iω be modifications made to s i
in the RT market in scenario ω ∈ Ω. Let C ap iω be the capacity of the bidder in P scenario ω ∈ Ω and C iω be the costs of generation. Denote S = i ∈N s i ,U ω = P P |Ω| |Ω| ω ω ω 1 2 1 2 i ∈N u i ,V = i ∈N v i , and u i = (u i , u i , . . . , u i ), v i = (v i , v i , . . . , v i ). Finally let p denote the price in the DA market and p ω denote the price in the RT mar-
ket for scenario ω. Then in a Nash-Cournot game, the i th generator’s profit maximization problem in such a two-settlement system can be formulated as follows: £ ¤ GenCoi (z −i )maximize p(S)s i + IEω p ω (U ω )u iω − peω (V ω )v iω −C iω (s i + u iω − v iω ) s i ,u i ,v i
subject to
s i + u iω − v iω ≤ C ap iω , u iω , v iω ≥ 0
60
∀ ω ∈ Ω,
where z i = (s i , u i , v i ). Finally our goal is to obtain a tuple {s i∗ , u i∗ , v i∗ }i ∈N such that ∗ ∗ ∗ s i∗ , u i∗ , v i∗ ∈ SOL(GenCoi (s −i , u −i , v −i )
The newsvendor game We consider here, a competitive extension of the canonical newsvendor problem in stochastic programming. Consider a noncooperative game between N newsvendors where the i th newsvendor would like to determine the number of newspapers to order, x i , in the face of uncertain demand. The newspapers are bought at cost c i (x i ) by the i th newsvendor through a private contract and sold P at a price p(Y ) where Y = i ∈N y i , and y i is the quantity sold by the newsvendor. We assume that the market is segmented and that every newsvendor has an access to a private uncertain demand d iω that he seeks to meet, where ω ∈ Ω.
Let y iω be the number of newspapers sold by the newsvendor in scenario ω. The
unsold newspapers may be returned to the manufacturer at a cost r iω (x − y iω ),
which depends on the realization of uncertainty. The i th newsvendor’s problem, parametrized by the other newsvendors’ orders is given by NVi (x −i )
£ ¤ maximize −c i (x i ) + IE p(Y )y iω + r iω (x − y iω ) x,y
y iω ≤ d iω
subject to
yω − x ≤ 0
x i , y iω ≥ 0
∀ ω ∈ Ω.
2.2.2 Outline of Literature In this section we outline some of the literature that precedes our work, focusing specifically on the Nash-Cournot model and its associated solution concepts. Nash-Cournot framework Nash equilibrium problems were born out of Nash’s ideas in 1950 [66] that clarified the distinction between cooperative and noncooperative games. von Neumann and Morgenstern [67] had provided the concept of a minmax solution to every game of finite strategies. Nash noted in [66] that this solution concept is permissible only if the players are allowed to communicate with each other and
61
arrive at a consensus in their bargaining. In the absence of any communication between players or any such arbitrating mechanism, the most reasonable solution one can assign to a game is that of an equilibrium solution; a tuple of the strategies of players from which no player has any incentive to deviate. Such a tuple of strategies is known as the Nash equilibrium for the noncooperative game, and has become the most widely applied solution concept in the domain of noncooperative games. Formally, if K i is the set of permissible strategies for player i ∈ N = {1, 2, . . . , n} and f i (x i ; x −i ) is the value to the player for playing strategy x i given strategies x −i played by other players, the Nash equilibrium is defined as x ∗ = (x 1∗ , x 2∗ , . . . , x n∗ ) such that ∗ ∗ f i (x i∗ ; x −i ) ≥ f i (x i ; x −i ) ∀ x i ∈ K i , ∀ i ∈ N.
Nash showed the existence of such equilibrium strategies in the domain of mixed strategies [66]. To understand what this means, consider an n−player game where K i is a finite set. A mixed strategy for player is a vector p i ∈
R|K | i
that satisfies pi ≥ 0
and
1T p i = 1.
A mixed strategy may be thought of as a probability distribution of a player’s preference over the elements of K i . The payoff of a mixed strategy p i is the ex£ ¤ pected payoff IE−i f i (p i ; x −i ) where the expectation is taken over p −i . A mixed strategy Nash equilibrium is hence a collection of randomization profiles, one profile for each player, from which no player has any incentive to deviate. Since then the Nash equilibrium and its refinements have come to be accepted as the most reasonable solution concept to any game that emerges from a noncooperative interaction between selfish agents. Tremendous research contributions have been made toward characterizing the existence and uniqueness of such equilibria in games with different kinds of mathematical properties. We refer the reader to the notes in section 1.9 of [27] for a review. An economic model of production and distribution where in the price (or the inverse demand function) of an item is a function of the total industry production of the item was first studied by Cournot in 1838 [20], and the models are often referred to as Cournot models. They form one of the most important production/distribution models [70]. Suppose we let N = {1, 2, . . . , N } be a set of producing firms and the quantity produced by firm i is denoted by x i . Let
62
X=
P
i ∈N x i .
The firm incurs costs c i (x i ) for producing x i units of goods. The
average market price p(X ) is dependent on the total quantum of goods produced, i.e. the model is a Cournot model. The firm i ’s problem is to determine x i∗ such that x i∗ ∈ arg max πi (x i ; x −i ) where πi (x i ; x −i ) := x i p(X ) − c i (x i ). Clearly such an x i∗ is dependent on the values x −i , the quantities produced by the other firms. The collection of n such problems can be regarded as a noncooperative game. A Nash equilibrium for such a game [90] is a tuple of strategies x ∗ = {x 1∗ , x 2∗ , . . . , x n∗ } such that x i∗ ∈ arg max x i p(x i +
X j 6=i
x ∗j ) − c i (x i )
∀ i ∈ N.
Such an equilibrium point has also been referred to as Nash-Cournot equilibrium or simply the Cournot equilibrium. Historically, the first line of research pursued of was that of theoretically proving the existence and the uniqueness of a Cournot equilibrium. The existence of such an equilibrium was first shown by Frank and Quandt [48] under the crucial assumptions that ∃ M > 0, such that p(X ) = 0 ∀ X ≥ M , concavity of πi and increasing cost functions c i . The first assumption implies that there is a finite quantity that the industry may produce so that the good is essentially free (‘the free good condition’). Under conditions stronger than that in [48], a new proof of existence and uniqueness of a Nash-Cournot equilibrium was provided in [90]. More recently in 1985, it was shown in [69] that if the free good condition holds and if each for i firm i ’s marginal revenue declines as the P aggregate j 6=i x j increases, then a Nash-Cournot equilibrium exists. Such a condition is implied by the concavity of p(X ). Finally in 1987 necessary and sufficient conditions for the uniqueness of the Nash-Cournot equilibrium were proved in [51]: if the Jacobian of the mapping π(x) = (π1 (x 1 ; x −1 ), π2 (x 2 ; x −2 ), . . . , πn (x n ; x −n )) evaluated at an equilibrium x ∗ is positive for any equilibrium x ∗ , then the equi63
librium is unique. The other line of research that followed these results is of computing NashCournot equilibria. Kolstand et al. [51] leveraged results from degree theory and linear complementarity theory to prove existence, and in a later paper [52] the same authors addressed the question of computing Nash equilibria using sequential linearization and Lemke’s algorithm. Spence [87] and Slade [86] provided some early conditions under which a Nash equilibrium problem, that of maximization of objectives of several firms, may be considered observationally equivalent to a mathematical program with a single objective function. Hence computing equilibria was shown to be equivalent to computing the solution of a mathematical program. Sherali et al. [63] developed a mathematical programming formulation in which the Nash-Cournot equilibrium problem was shown to be equivalent to a convex mathematical program. Since then the use of VIs, CPs and mathematical programs with equilibrium constraints (MPEC) to compute and characterize equilibria has been common. A few years after Sherali et al., Harker [36] provided a VI approach to the Nash-Cournot problem. Stackelberg-Nash-Cournot equilibria were later characterized by Sherali [85]. This chain of research continues to the present day. More references may be found in section 1.9 of [27]. Special attention may be given to the copious amounts of literature on using VIs and CPs particularly in the context of bidding in electricity markets. Due to the inherent uncertainty in the problem and Cournot nature of the competition, these problems result as stochastic Nash-Cournot equilibrium problems, like the ones described in section 2.2.1. Hobbs formulated the Nash-Cournot competition in electricity markets as a linear complementarity problem [41]. This approach was further extended in [62, 42]. Another MPEC approach to power markets can be found in [43]. More sophisticated models assuming network uncertainty and market power were discussed and analyzed by Oren et al. in [49]. Their formulations resulted in mixed linear complementarity problems for a single stage problem. However when the second stage problem is taken into account, they required the strategies to be subgame perfect, resulting in nonconvex MPECs for each agent, making the equilibrium problem an equilibrium program with equilibrium constraints (EPEC). In [97], the authors present an algorithm for solving such an EPEC. Indeed it is these game–theoretic subtleties in describing the underlying stochastic game that have been seen to severely impact the tractability of the models of oligopolistic competition in 64
electricity markets. The next section is devoted to a better understanding of these game–theoretic properties. Game–theoretic subtleties in stochastic games This section overviews three different notions of equilibria that emerge naturally in a stochastic game: Markov perfect equilibria, Oblivious equilibrium, S-adapted open loop equilibria. To introduce these we require the definition of an N -person discrete time infinite2 dynamic game. Definition 2.2.1 (N -person discrete time infinite dynamic game [6]) An N -person discrete time infinite dynamic game with prespecified fixed duration involves: 1. N = {1, . . . , N } is the set of players, 2. K = {1, . . . , K } is the set of stages of the game, where K is the maximum possible moves a player is allowed to make in the game. We allow K to be infinite. 3. An infinite set X with some topology defined on it, called the state set. The state of the game at stage k, x k ∈ X for all k ∈ K and for k = K + 1, x k+1 is called the final state of the game. 4. A collection of infinite sets Uki with some topology defined on them for each i ∈ N and each k ∈ K. The set Uki is the set of actions at stage k for player i . 5. A function F k : X ×Uk1 × . . . ×UkN → X defined for each k ∈ K and satisfies x k+1 = F k (x k , u k1 , . . . , u kN ) ∀ k ∈ K The above equation is the state equation. Note that this equation captures the memory of the actions of the players in the stages preceeding the current stage implicitly through the current state x k . 6. A set Yki with a topology defined on it defined for each k ∈ K and each i ∈ N called the observation set of player i at stage k. 2
The game is called an infinite dynamic game since the strategy sets are assumed to be infinite.
65
7. A function h ki : X → Yki defined for each k ∈ K and each i ∈ N is the measurement function for player i at stage k. The observed state of the game for player i at stage k, y ki is given as y ki = h ki (x k )
k ∈ K,
i ∈ N.
8. A finite set I ki , defined for each k ∈ K and i ∈ N as a subcollection of 1 2 N {y 11 , . . . , y k1 ; y 12 , . . . , y k2 ; . . . ; y 1N , . . . , y kN ; u 11 , . . . , u k−1 ; u 12 , . . . , u k−1 ; . . . ; u 1N , . . . , u k−1 }
determines the information gathered and recalled by player i at stage k, is called the information set. The sequence of N −tuples {(I k1 , . . . , I kN )} is the information pattern of the game. 9. Each information set resides in an information space, denote by I ki which is an appropriate subset of ((Y11 × . . . × Yk1 ) × (Y12 × . . . × Yk2 ) . . . (y 1N × . . . × y kN )× 1 2 N (U11 × . . . ×Uk−1 ) × (U12 , . . . ,Uk−1 ) . . . (U1N , . . . ,Uk−1 )).
(2.1)
10. Mappings γik : I ki → Uki belonging to a class Γik defined for each i ∈ N and k ∈ K are the set of permissible strategies. Note that strategies differ from actions. Strategies are maps that take information to action. γik is refered to as the strategy for player i at stage k. The tuple γi := (γi1 , . . . , γiK ) is the strategy for player i for the game. The space of all mappings γi is the strategy space of player i 11. Finally a functional L i : (X ×U11 × . . .U1N ) × (X ×U21 × . . .U2N ) × . . . (X ×UK1 × . . .UKN ) → ℜ is the cost (objective function) of player i .
Most elements in the above definition are self-explanatory and may be likened to those found in static games. Special attention should be paid to the definition of strategies in a dynamic game. A dynamic game maintains a clear distinction between a strategy and an action. A point in the action set is an action. A
66
strategy is a mapping from the information set to the set of actions. It is hence an extensive plan of actions based on the information available to the player. In a dynamic game, information is gathered by a player through the evolution of the stages of the game, which influences the strategies of the player. A player’s strategy is hence most explicitly described as a function on the information set. Evidently, the dynamic game (and it’s special case, the stochastic game) leads to rich classes of equilibria. The following sections elucidate some of these concepts. 1. Markov perfect equilibrium [60, 61] : With a minor modification in the above definition, if K is taken to be ∞, one obtains the infinite horizon game. The Markov perfect equilibrium concept introduced by Maskin and Tirole [60, 61] is relevant in the context such infinite horizon games to analyze the long run behaviour of firms in competition. The model (for two players) is as follows. The game is played sequentially; player 1 moves at odd instants of time and player 2 moves at even instants of time. Since the long run behaviour is in question, which firm plays first is immaterial. The action sets of players Uki are assumed to be bounded and for simplicity assumed to be the same for all stages k. At any instant t , the profit of each firm is assumed to be a function the previous actions, but not of the time. i.e. πi = πi (u t1 , u t2 ) There is a common discount factor associated with each firm δ = exp(−r T ) where is T is the time between two periods t i , t i +1 . Each firm’s intertemporal profit is assumed to be additive, Li =
∞ X
δπi (a t1 , a t2 ).
t =1
Markov perfect strategies are required to be representative of the long term behaviour of firms. Hence the natural requirement for such strategies is that they should not depend on the entire history of actions. The requirement imposed by Maskin and Tirole is that of ‘Markovian’ behaviour: the strategies depend only on the most recent action of the opponent. In light of the definition of a dynamic game the information sets consist of only the most recent strategy of the opponent. Hence a Markovian strat-
67
egy is a mapping γi : U −i → U i . The resulting actions are 2 1 1 2 u 2k = γ2 (u 2k−1 ), u 2k+1 = γ1 (u 2k ).
Strategies γ1 and γ2 are said to form a Markov perfect equilibrium if and only if 2 1 1 (a) u 2k = γ2 (u 2k−1 ) maximizes L 2 given action u 2k−1 and assuming that
each firm plays according to γi in the future. 1 2 2 (b) u 2k+1 = γ1 (u 2k ) maximizes L 1 given action u 2k and assuming that
each firm plays according to γi in the future. In the context of a stochastic game, the summation in the definition of L i is replaced by the expectation over the time horizon, while strategies γi are random functions and each action is a random variable. 2. Oblivious Equilibrium [95]: The oblivious equilibrium concept is an alternative to the Markov perfect equilibrium concept for analyzing stochastic dynamic games. In [95], the authors define the additional quantity, the state of a player, which determines his ability to compete in the game. The state of a player i at time t is denoted by x i ,t ∈ N. The set of actions is assumed to be common for all players and the same at each time instant, denoted as U . The state of the game is a vector over individual states specifying the number of players in each state. x −i ,t denotes the state of his competitors. A strategy is a mapping γ : N×Nn−1 → U . The space of all such strategies Γ is supposed to be shared by all players. Thus the action of player i is u i ,t := γ(x i ,t , x −i ,t ). The state of the player transitions based on a random process ζi ,t , his action and a common function F x i ,t +1 = x i ,t + F (u i ,t , ζi ,t +1 ). An agent’s single period profit is denoted by π(x i ,t ; x −i ,t ). L i is assumed to have an additive form, with discount factor β. The value function is taken to be the net present value for an agent at state x when it’s competitors are in state s, given that it’s opponents all use strategy γ while it uses strategy γ0 is
" 0
V (x, s|γ , γ) = IE
∞ X
#
β
k−t
0
π(x i ,k , s −i ,k )|x i ,t = γ , x −i ,t = γ
k=t
68
The kind of equilibrium sought is symmetric, such that all players use a common strategy. An oblivious equilibrium comprises a strategy γ ∈ Γ such that sup V (x, s|γ0 , γ) = V (x, s|γ) γ0 ∈Γ
∀x ∈ N, x −i ∈ Nn−1 .
3. S-adapted open loop equilibria[37]: It may be noted that the equilibrium of the stochastic Nash-Cournot game in definition 2.1.1 is neither a Markov perfect equilibrium nor an oblivious equilibrium. The game in 2.1.1 is a peculiar case of the dynamic game 2.2.1. The stochastic NashCournot game may be viewed as a dynamic game with two stages and two sets of decisions (x i , y i ) for each player. For a finite scenario set Ω the game in it’s extensive form can be represented as a tree. Note however that the decision y iω taken at the ωth node of the second stage, is taken with the same information about the game as the decision x i . i.e. the game, though with two stages, has indeed only one information set, denoted by I . A strategy for player i is hence a mapping γi = (x i , y i1 , y i2 , . . . , y i|Ω| ) : I → K i ×
Y
K iω (x i ),
such that x i ∈ K i , y iω ∈ K iω (x i ). Hence the stochastic Nash-Cournot game despite having two stages, is static in nature, since the information sets of the players are static. Since the strategies y iω are determined without any extra information about the first stage equilibrium, they cannot be expected to be subgame perfect. Noting this, we now come to a third type of equilibrium that arises by treating the stochastic Nash-Cournot game as a special case of a stochastic dynamic game. We source the following definition from [6]. Definition 2.2.2 (Closed loop strategy) Consider a two stage dynamic game with action sets U1i ,U2i where U2i = U2i (1) ×U2i (2) × . . . ×U2i (|Ω|) for some finite set Ω. Let K iω (u) ⊆ U2i (ω) for each u ∈ U1i and each ω ∈ Ω be given sets. A closed loop strategy for player i for such a game is a pair (γi1 , γi2 ) such that ∃u 1i ∈ U1i such that γi1 (z) = u 1i for all z ∈ I 1i and γi2 = (γi2,1 , γi2,2 , . . . , γi2,|Ω| ) is
69
a vector of mappings such that γi2,k : U1i → U i (k) s.t . γis,k (u) ∈ K iω (u) ∀ u ∈ U1i . In effect definition 2.2.2 says that a closed loop strategy is a pair of mappings. The first mapping is a constant function. The second is a finite collection of functions that map the first stage action set to the second stage action set. This means that the information set of the second stage is I 1i together with the action of the first stage u 1i . In addition each function in the collection is required to be feasible with respect to a set parametrized by the first stage action. An open loop strategy in comparison, is a pair (γi1 , γi2 ) such that mappings γi2 are completely independent of U1i . The reader may verify that if the SNC game were to be presented as two stage dynamic game then the strategies (x i , y i ) would form a closed loop pair. Indeed, the SNC game is not a dynamic game in the strictest sense, since the second stage decisions y i are taken simultaneously with x i , with no new information about the x i . The strategies (x i , y i ) do not form an open loop pair either, since y iω is required to be feasible with respect to K iω (x i ). Haurie et al. [37] show that due to the feasibility requirement y iω ∈ K i (x i ) in the SNC game, it can be viewed as a special kind of dynamic game. They define the following strategy type that is intermediate between closed loop and open loop strategy. Definition 2.2.3 (S-adapted open loop strategy [37]) Consider a two stage stochastic dynamic game as in definition 2.2.2. A strategy pair (γi1 , γi2 ) for player i is called S-adapted open loop strategy if ∃ u 1i ∈ U1i such that γi1 (z) = i i u 1i for each z ∈ I 1i and ∃ u 2,ω ∈ K iω (u 1i ) ∀ ω ∈ Ω such that γi2,ω (u) = u 2,ω for i each u ∈ U1i such that u 2,ω ∈ K iω (u), for each ω ∈ Ω.
An S-adapted open loop strategy is thus a pair of constant mappings, such that the second mapping is a collection that is feasible with respect to the first stage mapped value. It is the constancy of the feasible second stage mappings that makes them independent of the first stage mappings. It is easy to see the any feasible strategy for the stochastic game is an S-adapted open loop strategy. Since the equilibrium sought by the stochastic Nash-Cournot game is in S-adapted open loop strategies, the equilibrium will be called the S-adapted open loop equilibrium. 70
2.3 The Deterministic Nash-Cournot Problem We consider problems at the interstices of stochastic programming and gametheory. In the instance of the former, planners are faced with making decisions in the face of an uncertain future. Game-theoretic problems on the other hand imply that mere optimizing will not suffice. Such problems can be solved as variational inequalities or complementarity problems [27]. However, when the agent problems are two-period stochastic programs, then the resulting equilibrium problems are stochastic variational inequalities or complementarity problems. In this section we formulate the deterministic Nash-Cournot problem and prove some properties of this problem. In the next section, we leverage properties of the deterministic problem to formulate its stochastic version. Throughout this section, N = {1, 2, . . . , N } denotes the set of players, K i is the set of feasible decisions of player i and Ω (a finite set) denotes the sample space of the Q random variable ω. We use K = K i . Definition 2.3.1 (DNC) The deterministic Nash-Cournot game is a single period game in which the i th player solves the deterministic convex program given by Agent(x −i )
maximize p(X )x i − f i (x i ) xi
subject to
c i (x i ) ≥ 0 x i ≥ 0,
where f i and c i are convex and concave twice continuously differentiable functions. A Nash-Cournot equilibrium of the above game is given by a tuple x ∗ = {x 1∗ , x 2∗ , . . . , x N } where ∗ x i∗ ∈ arg max Agent(x −i )
Note that in (DNC), x i ∈
∀i ∈N
R for each i . The price function is expected to satisfy
the following collection of assumptions: Assumption 2.3.2 The price function p(X ) satisfies the following requirements: 1. The function p(X ) is a non-increasing twice-continuously differentiable function. 2. The function X p(X ) is a strictly concave function of X ≥ 0. 71
Cost and constraint functions are assumed to satisfy the following: Assumption 2.3.3 The cost and constraint functions f i (x i ) and c i (x i ) are twice continuously differentiable functions for all i ∈ N. Furthermore, there exists a positive constant κ such that f i00 (x) > κ for all x ∈ K i and for all i ∈ N. Lemma 2.3.4 Suppose X p(X ) is a strictly concave function in X . Then if p(X ) is either non-increasing or convex, then the the concavity of x i p(X ) for all i follows. Proof : We reproduce this proof from [63]. Note that x i p(X ) = x i p(x i + Assume that the required result does not hold and r 00 (x i ) = (x i p(x i + X
P
j 6=i −i 00
x i ).
)) > 0.
Then we have r 00 (x i ) = 2p 0 (x i + X −i ) + x i p 00 (x i + X −i ) = 2p 0 (x i + X −i ) + (X −i + x i )p 00 (x i + X −i ) − X −i p 00 (x i + X −i ) < −X −i p 00 (x i + X −i ). If r 00 (x i ) is not strictly concave, it follows that for nonnegative X −i , p 00 (x i + X −i ) < 0. However, if p(x i + X −i ) is a convex function, this results in a contradiction. Alternately, if we assume that p(x i + X −i )0 ≤ 0, then it follows that x i p 00 (x i + X −i ) > 0, implying that p 00 (x i + X −i ) ≥ 0. But this contradicts, p 00 (x i + X −i ) ≤ 0. It follows that r (x i )00 ≤ 0. It follows that each agent problem is essentially a problem with a strictly concave objective and constraints defined a closed and convex set. The optimality conditions of this agent problem are given by F i (x i )T (y i − x i ) ≥ 0, where
∀y i ∈ K i ,
F i (x i ) = f 0 (x i ) − p 0 (X )x i − p(X ), K i = {y i : c i (y i ) ≥ 0, y i ≥ 0}.
The resulting VI corresponding to the game-theoretic problem is given by F(x)T (y − x) ≥ 0,
72
∀y ∈ K ,
where F(x) is given by F(x) = G(x) + R(x), −p 0 (X )x 1 − p(X ) f 10 N Y . .. and K := .. , R(x) := Ki . G(x) := . i =1 −p 0 (X )x N − p(X ) f N0
The Jacobian of the mapping F is given by ∇F = ∇G+∇R. The strict monotonicity of G follows immediately from the strict convexity of the objectives of each agent. In addition, the Jacobian of R is defined as
HR .11 . ∇R = . R HN 1
R H1N .. . R . . . HN N
... .. .
H Rjj = −2∇ j p(X ) − x j ∇2j j p(X ) = H Rj` − ∇ j p(X ) H Rj` = −∇ j p(X ) − x j ∇2j ` p(X ) H`Rj = −∇ j p(X ) − x ` ∇2` j p(X ). Specifically, we can write the second-derivatives of R as H Rjj = β + α j
H Rj` = α j
H`Rj = α` ,
where β = −p 0 (X ) and α j = −p 0 (X ) − x j p 00 (X ). This allows us to show to prove the following existence and uniqueness result for the associated variational inequality. Note that this result does not require the concavity of of p(X ) or the free good assumption of [69], cited in section 2.2.2. Theorem 2.3.5 Suppose the price, cost and constraint functions in (DNC) satisfy assumptions 2.3.2 and 2.3.3. Then the Nash-Cournot equilibrium problem has at most one solution. Proof : Given the strict convexity of f i and strict concavity of c i , the matrix ∇G is clearly positive definite implying that ∇G is a P-matrix. To show that ∇R is a P matrix, we provide a proof based on determinants
73
sourced from [53]. Recall that ¯ ¯β + α α1 α1 ¯ 1 ¯ ¯ α2 β + α2 α2 ¯ det H N = ¯ . .. .. ¯ .. . . ¯ ¯ ¯ αN αN αN
¯ ¯ ¯ ¯ ... α2 ¯¯ .. ¯¯ .. . . ¯ ¯ . . . β + αN ¯
α1
...
It suffices to show that all principal minors of this matrix are positive. By adding all the rows to the first one and subtracting from the first column, we have that ¯ ¯β + P N α ¯ i =1 i ¯ ¯ α2 ¯ det H N = ¯ .. ¯ . ¯ ¯ ¯ αN
= βN −1 (β +
N X
¯ 0 ¯¯ ¯ β 0 . . . 0 ¯¯ .. .. . . .¯ . .. ¯¯ . . ¯ 0 0 . . . β¯
0 0 ...
αi )
k=1
= (−p 0 (X ))N −1 (−(N + 1)p 0 (X ) +
N X
−x k p 00 (X )).
k=1
Recall that (r i (x i ; x −i )00 = 2p 0 (X ) + x i p 00 (X ) implying that N X
(−p 0 (X ))N −1 (−(N + 1)p 0 (X ) +
−x k p 00 (X ))
k=1 N X
= (−p 0 (X ))N −1 (−(N − 1)p 0 (X ) +
−r k (x k ; x −k )00 )
k=1
> 0, where p 0 (X ) < 0 and r j (x j ; x − j )00 < 0, j = 1, . . . , N . It follows that det H N > 0. It can be further shown that every principal minor is positive. It follows that the ∇F is a P-matrix implying that F is a P-function. Such a property ensures that VI(K,F) has at most one solution. Proposition 2.3.6 Suppose the price and cost functions satisfy assumptions 2.3.2 and 2.3.3. Then an equilibrium to the Nash-Cournot equilibrium problem exists and is unique. Proof : Observe that K =
Q
K i . It is easy to see that K is closed. Two cases arise
74
(a.) K is bounded. It is hence compact. Then by Kakutani’s fixed point theorem there exists a solution to the VI(K , F). (b.) K is not bounded. i.e. ∃i¯ ∈ {1, . . . , N } such that K i¯ is unbounded. WLOG we assume i¯ = 1. For case (b.) our proof relies on noting that our variational problem can be cast as a VI on a Cartesian product of closed convex sets: y 1 − x1 0 F 1 (x)T . . . ≥ .. , ∀y i ∈ K i , .. .. T yN − xN 0 F N (x)
i = 1, . . . , N .
Given such a characterization, a weaker sufficient condition leads to existence of an equilibrium. Specifically, we need to show that that there exists an x ref ∈ K such that for some j ∈ {1, . . . , N } we have lim
x j ∈K j ,kx j k→∞
F jT (x j ; x − j )(x j − x ref j )≥0
We show that this holds for j = i¯ = 1. We begin by noting that F 1 (x) = f 10 (x 1 ) − p 0 (X )x 1 − p(X ). Take some x ref ∈ K . Denote δx = X − x 1 , δref = X ref − x 1ref lim
x 1 ∈K 1 ,kx 1 k→∞
= =
F 1 (x)T (x 1 − x 1ref )
lim
(F 1 (x) − F 1 (x ref ) + F 1 (x ref ))T (x 1 − x 1ref )
lim
F 1 (x ref )(x 1 − x 1ref ) +
x 1 ∈K 1 ,kx 1 k→∞
x 1 ∈K 1 ,kx 1 k→∞
lim
[( f 10 (x 1 ) − f 10 (x 1ref )(x 1 − x 1ref )
x 1 ∈K 1 ,kx 1 k→∞
− (p(X ) − p(X ref ))(x 1 − x 1ref ) − (x 1 p 0 (X ) − x 1ref p 0 (X ref ))(x 1 − x 1ref )]. Note also that x 1 p 0 (x 1 + δx ) − x 1ref p 0 (x 1ref + δref ) = x 1 p 0 (x 1 + δx ) − x 1ref p 0 (x 1ref + δx ) + x 1ref p 0 (x 1ref + δx ) − x 1ref p 0 (x 1ref + δref ). By mean value theorem and assumptions 2.3.3, 2.3.2, ∃x˜1 , x˜2 , x˜3 and c˜1 > 0, c˜2 ≥ 0, c˜3 ≥ 0,
75
such that f 10 (x 1 ) − f 10 (x 1ref ) = f 100 (x˜1 )(x 1 − x 1ref ) > c˜1 (x 1 − x 1ref ), −(p(X ) − p(X ref )) = −p 0 (x˜2 )(x 1 + δx − x 1ref − δref ) ≥ c˜2 (x 1 − x 1ref + δx − δref ), −(x 1 p 0 (x 1 + δx ) − x 1ref p 0 (x 1ref + δx )) = −(xp 0 (X ))0 (x˜3 )(x 1 − x 1ref ) ≥ c˜3 (x 1 − x 1ref ). It is clear that F 1 (x)T (x 1 − x 1ref ) > c(x 1 − x 1ref )2 + O((x 1 − x 1ref )) + O(1), where c = c˜1 + c˜2 + c˜3 . Hence lim
x 1 ∈K 1 ,kx 1 k→∞
F 1 (x)T (x 1 − x 1ref ) = ∞
and the existence result follows. Uniqueness follows from theorem 2.3
2.4 The Stochastic Nash-Cournot Equilibrium Problem In this section, we provide a formal description of the stochastic Nash-Cournot problem. Such a problem is set in a market which has two settlements. The first settlement is a consequence of a game between sellers to determine the price of electricity, with a pre-selected price function. The second settlement occurs just before the physical trading of the asset and represents recourse actions taken by agents, under a variety of realizations of uncertainty. The stochastic Nash-Cournot equilibrium problem of interest is a static game in which agents make 1. First-period decision that specifies the quantities (and therefore the price) committed by each agent and 2. Second-period recourse decisions that represent recourse actions taken by agents in accordance with high fuel prices or depleted capacity. In the single-period game, this results in agents maximizing the sum of firststage profits and expected second-stage profits, arising from recourse actions.
76
Expectedly, a multitude of formulations may arise given such a general framework. Our intent in this section is to present one formulation and provide some direction on where modifications may arise.
2.4.1 Formulating the Stochastic Nash-Cournot Problem As mentioned in the earlier section, the stochastic Nash-Cournot problem requires each player to maximize his expected profit over a measure IP. This measure is specified over the randomness in the second period and can take two forms: 1. Common information: The distribution IP specifies the likelihood that a specific commodity is going to be priced in the future. For instance, the random variable specifying the price of oil could be one such commodity, the assumption being that the distribution of oil prices is known to all the players. 2. Private information: In many instances, players may have private information that is not shared across players. For instance, such information could pertain to the likelihood of certain production units being low on capacity. Often such information has to be made available to the market authority. In this chapter we work with only a common information structure among agents. We suggest a Nash-Cournot framework that offers significant breadth in terms of the type of games played in both periods. While the first-period game is a game in quantities with a well-defined price function, the second-period game is a game parametrized by first period decisions. This problem is also taken to be a Cournot game in which agents determine deviation levels from the prescribed Cournot bids in the first period. Consider an equilibrium problem in which players are competing in a twoperiod market in which the prices are specified by p(X ) and p ω (U ω ) where X=
N X
xi ,
Uω =
i =1
N X i =1
u iω .
We assume that all the price functions satisfy assumption 2.3.2.
77
Definition 2.4.1 (SNC) Consider a game in which the i th player solves the parametrized stochastic optimization problem given by
SNC-Ag(x −i , u −i )
maximize
e ω )v iω [p(X )x i − f i (x i )] + IEω [p ω (U ω )u iω − p(V
−h iω (u iω ) + r iω (v iω )]
x i ,u i ,v i
subject to
c i (x) ≥ 0 ω ω ω σi (x i ) + s i (u i ) ≥ 0 τω (x i ) + t iω (v iω ) ≥ 0 i x i , u iω ≥ 0,
∀ω ∈ Ω.
where f i , c i , h iω , s iω , t iω , τω and σω represent twice continuously differentiable funci i tions. The stochastic Nash-Cournot equilibrium is given by {x i∗ , u i∗ }i ∈N where ∗ ∗ (x i∗ , u i∗ ) ∈ SOL(SNC-Ag(x −i , u −i )),
∀ i ∈ N.
Recall that such an equilibrium falls under the paradigm of S-adapted open loop equilibria of stochastic games as shown in section 2.2.1. We conclude this subsection with existence statements for the equilibria of (SNC). Proposition 2.4.2 Consider the game (SNC). Assume that the cost and constraint functions satisfy assumption 2.3.3. Specifically, {h iω }ω∈Ω,i ∈N and { f i }i ∈N are con-
vex; the constraint functions c i and {τω , t ω , σω , s iω }ω∈Ω,i ∈N are all concave. Supi i i
pose the price functions p(·), p ω (·) satisfy 2.3.2, and peω (V ω )v iω is convex in v iω . Then a stochastic Nash-Cournot equilibrium exists. Proof : Consider the problem (SNC). If ω is assumed to be defined on a discrete probability measure of size |Ω|, then the agent problems are merely larger con-
vex programs. Specifically, the i th agent is maximizing over the set of variables:
Ã
xi ui
!
v i1 u i1 . . . . , ui = . , vi = . u i|Ω| v i|Ω|
78
with cost and revenue functions given by f¯i (x i , u i ) :=
f i (x i ) +
|Ω| X j =1
and ρ¯ i (x i , u i , v i ) := p(X )x i +
j
j
j
j
$ j (h i (u i ) − r i (v i ))
Ω X j =1
$ j (p j (U j )u i − pe j (V ω )v i ), j
j
where $ j = Prob{ω = j }, j ∈ Ω. The convexity of the cost function f¯ and the concavity of the revenue function ρ¯ follows from assumption 2.3.2 and the hypothesis of the theorem. Each agent is optimizing over a closed convex set K i where ω ω K i = {(x i , u i ) : c i (x i ) ≥ 0, σω i (x i ) + s i (u i ) ≥ 0, ω ∈ Ω}.
The resulting VI(F, K ) is specified over the Cartesian product K =
Q
i ∈N K i
and
the mapping F equates to ∇z1 f¯1 (z 1 ) − ∇z1 ρ¯ 1 (z 1 ; z −1 ) ³ ´ .. , zi = xi , ui , v i . F(z) = . −N ¯ ∇z f N (z N ) − ∇z ρ¯ N (z N ; z )
N
N
Fi (z) is given by
f i0 (x i ) − p 0 (X )x i − p(X )
$1 (h i10 (u i1 ) − p 10 (U 1 )u i1 − p 1 (U 1 )) .. . |Ω|0 |Ω| |Ω| |Ω| |Ω|0 |Ω| |Ω| |Ω| Fi (z) = $ (h i (u i ) − p (U )u i − p (U )) . $1 (−r i10 (v i1 ) + pe10 (V 1 )u i1 + pe1 (V 1 )) .. . $|Ω| (−r i|Ω|0 (v i|Ω| ) + pe|Ω|0 (V |Ω| )v i|Ω| + pe|Ω| (V |Ω| ))
We see that Fi has the same properties as the corresponding mapping in (DNC). Once again if K is compact, then existence follows from Kakutani’s fixed point theorem. If K is unbounded, it can be inferred from the proof of proposition 2.3.6 that ∃ x ref ∈ K such that lim
x i ∈K i ,kx i k→∞
F i (x)T (x i − x 1ref ) > 0 for somei ∈ N
79
To prove uniqueness, we rearrange the vector z in the following manner: ³ ´T |Ω| 1 |Ω| 1 1 z = x 1 , . . . , x N , u 11 , . . . , u N , , . . . , u 1|Ω| , . . . , u N , v1 , . . . , v N , . . . , v 1|Ω| , . . . , v N
and get F as follows.
f 10 (x 1 ) − p 0 (X )x 1 − p(X ) .. .
f N0 (x N ) − p 0 (X )x N − p(X ) 1 10 1 10 1 1 1 1 $ (h 1 (u 1 ) − p (U )u 1 − p (U )) .. . 1 10 1 10 1 1 1 1 $ (h N (u N ) − p (U )u N − p (U )) .. . |Ω| |Ω|0 |Ω| $ (h 1 (u 1 ) − p |Ω|0 (U |Ω| )u 1|Ω| − p |Ω| (U |Ω| )) .. F= . |Ω| |Ω|0 |Ω| |Ω| $ (h N (u N ) − p |Ω|0 (U |Ω| )u N − p |Ω| (U |Ω| )) $1 (−r 110 (v 11 ) + pe10 (V 1 )v 11 + pe1 (V 1 )) .. . 1 10 1 10 1 1 1 1 $ (−r N (v N ) + pe (V )v N − pe (V )) .. . |Ω| $ (−r |Ω|0 (v |Ω| ) + pe|Ω|0 (V |Ω| )v |Ω| + pe|Ω| (V |Ω| )) 1 1 1 .. . |Ω|0 |Ω| |Ω| |Ω|0 |Ω| |Ω| |Ω| |Ω| $ (−r N (v N ) + pe (V )v N − pe (V ))
The Jacobian of F is given by ∇F =
A B1 ..
. B |Ω| C1 ..
. C |Ω|
80
,
where all the offdiagonal blocks are 0. The submatrices A, B i ,C i , i ∈ Ω are defined as
α11
A=
αN 1
... .. .
α1N
. . . αN N
Bi = βiN 1
βi11
... .. .
βi1N
. . . βiN N
Ci =
γi11 γiN 1
... .. .
γi1N
. . . γiN N
where ∀ i , j ∈ N, i 6= j and ∀ k ∈ Ω, αi i = f i00 (x i ) − p 00 (X )x i − 2p 0 (X ) αi j = −p 00 (X )x i − p 0 (X ) α j i = −p 00 (X )x j − p 0 (X ) βkii = $k (h ik0i (u ik ) − p k00 (U k )u ik − 2p k0 (U k )) βkij = $k (−p k00 (U k )u ik − p k0 (U k )) βkji = $k (−p k00 (U k )u kj − p k0 (U k )) γkii = $k (−r ik0i (u ik ) + pek00 (U k )u ik + 2pek0 (U k )) γkij = $k (pek00 (U k )u ik + pek0 (U k )) γkji = $k (pek00 (U k )u kj + pek0 (U k )) Then by using lemma 2.3.4, we claim by the concavity of the revenue functions ρ i that the matrices A, B ω , ω ∈ Ω are each P-matrices. Let A be a principal submatrix of ∇F. Clearly A is a block diagonal matrix, comprising of blocks Ak , say. For each k, Ak is a principal submatrix of one of {A, B 1 , . . . , B |Ω| ,C 1 , . . . ,C |Ω| }. Hence det(Ak ) > 0 for all k since A, B i ,C i are P-matrices. By Cramer’s rule, Y det(A) = det(Ak ) > 0. So ∇F is a P-matrix. Therefore, the function F is a Pk
function implying that the resulting VI has at most one solution. This together with the proof of existence implies that VI(F, K ) has a unique solution.
81
3
An Inexact-Newton Splitting (INS) Method A pressing question that accompanies the construction of stochastic equilibrium problems is the development of algorithmic schemes. Importantly such schemes should be characterized by the following properties: 1. Globally convergent behavior: This is a property that guarantees that from any starting point, the algorithm can construct a sequence of iterates whose limit point is the required equilibrium point. 2. Rate of convergence results: Such rate specifications specify the rate of convergence of the iterates close to the solution. 3. Scalability: Stochastic equilibrium problems can grow with the size of the probability support. In keeping with the paradigm of optimization algorithms for stochastic programming problems, our intent is to propose a scheme whose effort grows linearly with the number of realizations that the random variable can take. In this section, we describe an algorithm that possesses the required properties of global convergence and scalability. Furthermore, we discuss the local convergence behavior of the overall scheme as well as that of the subproblem solver. In section 5.1, we outline the basic Newton framework we employ. Standard Newtonian schemes need to be globalized and we achieve this by relying on existing research that uses a D-gap merit function for precisely such a purpose. Unfortunately, in practical settings, such a scheme cannot be employed unless one can solve the Newton problem efficiently. Specifically, can one design a scheme such that the effort for solving the linearized problem grows linearly with the size of the probability measure? Section 5.2 is devoted to such a goal which we attain through the use of matrix splitting techniques. Convergence of the overall scheme and underlying subproblem is studied in section 5.3. In section 5.4, we provide some generalizations to the Nash framework that still allow for employing our method. This requires proving that the 82
resulting linearized problem satisfies certain properties. Finally, in section 5.5, some numerical evidence of the scalablility of the approach is given.
3.1 An Inexact Newton Framework Our basic problem of interest is the Nash equilibrium problem. Such a problem can be broadly captured by a setting where each agent solves a problem given by (SNC). We begin with a deterministic restriction of the problem. Consider a Nash equilibrium problem in which the i th agent solves a problem given by minxi f (x i ; x −i )
A(x −i )
s/t
c i (x i ) ≥ 0
(λi )
x i ≥ 0,
where f i and c i are convex and concave twice continuously differentiable functions, respectively. The resulting sufficient conditions are given by the complementarity problem
x1 ∇ f 1 − ∇c 1T λ1 . .. . . . x ∇ f − ∇c T λ N N N N 0≤ ⊥ ≥ 0. λ1 c 1 (x i ) . .. .. . λN c N (x N )
Note that in the case of the Nash-Cournot framework, f i (x i ; x i ) = p(X )x i −h i (x i ), where p(X ) is the price function and h i (x i ) is the cost of producing x i units. The aforementioned complementarity problem can be compactly written as 0 ≤ z ⊥ F(z) ≥ 0,
83
where F is a continuously differentiable multifunction defined as
∇ f 1 − ∇c 1T λ1 x1 . . . .. . ∇ f − ∇c T λ x N N N N F(x, λ) := and z = . λ1 c 1 (x i ) . . .. .. c N (x N ) λN
This problem can be equivalently written as the variational inequality: (y − z)T F(z) ≥ 0,
∀y ∈
Rn¯ .
We propose an inexact Newton method for solving such a problem. An exact Newton method requires solving a sequence of linearized problems, given by (y − z)T Fk (z) ≥ 0,
∀y ∈
Rn¯ ,
where Fk (z) = F(z k ) + ∇F(z k )(z − z k ). The matrix ∇F(z k ) can be written as ∇F(z k ) =
Hk
à Hk
−∇c T
!
∇c 0 P ∇211 f 1 − λ1i ∇2 c 1i = ∇21N f N
... .. .
∇21N f 1
. . . ∇2N N f N −
P
λN i ∇2 c N i
(3.1)
It follows that this variational inequality is equivalent to a complementarity problem: 0 ≤ z ⊥ ∇Fk z + Fk ≥ 0. The solution to this linearization is termed as z k+1/2 in keeping with the notation from [27]. However, such a solution provides a direction along which one attempts to obtain descent with respect to a suitable merit function. Given that our subproblems can be massive linear complementarity problems, we choose to solve them inexactly. This is operationalized through the
84
use of a sequence of nonnegative scalars {η k }. Moreover, an inexact solution of the complementarity problem would require an z k+1/2 such that Z k+1/2 (Fk + ∇Fk z k+1/2 ) ≤ η k e
(3.2)
where Z k+1/2 = diag(z k+1/2 ). Our merit function of choice is the D-gap merit function introduced in [71] and is defined as follows. Definition 3.1.1 (10.3.1. from [27]) Given a mapping F : closed convex set in
R
n¯
Rn¯ → Rn¯ and K
a
. Let a and b be scalars such that b > a > 0. Then the
D-gap function of VI(K , F) is defined as θab (z) ≡ θa (z) − θb (z), where
z∈
Rn¯ ,
´ ³ c θc (z) = sup F(z)T (z − y) − (z − y)T G(z − y) , 2 y∈K
and G is a symmetric positive definite matrix. The D-gap merit function is indeed an unconstrained merit function of VI(K , F) in that it is nonnegative over all the entire set K and θab (z) = 0 ⇔ z ∈ SOL(K , F). While the original Newton method (using exact solutions to the subproblem) was provided in [71], we adhere to a version provided in [27]. Importantly, the method suggested in the latter reference employs weaker assumptions in proving the global convergence of the scheme. Importantly, our scheme relies on inexact solutions and we provide a corresponding extension to the stated convergence results, under such an assumption. We now state the algorithm in its complete form: The key difference between our approach and that in [27] is the technique for accepting a direction. While in both references [71, 27], if the solution to the subproblem does not result in sufficient descent, then one immediately obtains a descent direction using the negative of the gradient of the merit function.
85
Algorithm 4: An inexact Newton scheme for nonlinear CPs 0 n¯ 0 initialize z ∈ , 0 < a < b, ρ > 0, r > 1, γ ∈ (0, 1), η¯ < 1, jˆ ∈ N and a
R
1 2
3
4
sequence of nonnegative scalars {η k }; k = 0; if z k ∈ SOL(K , F) then stop end Let z k+1/2 ∈ FEA(K , Fk ) and Z k+1/2 (Fk + ∇Fk z k+1/2 ) ≤ η k e and d k = z k+1/2 − z k ; Set j = 1; while ∇θab (z k )T d kj > ρkd kj kr or j ≤ jˆ do j
j −1
5
Let η k = η¯ η k ;
6
Let z k+1/2 ∈ FEA(K , Fk ) and (Z jk+1/2 )(Fk + ∇Fk z k+1/2 ) ≤ η k e, and j j
j
d kj = z k+1/2 − zk ; j 7
8
j = j + 1; end if ∇θab (z k )T d kj > ρkd k kr then d k := −∇θab (z k ) else d k = d kj .
9
10 11
end Find the smallest i k ∈ Z+ such that θab (z k + 2−i k d k ) ≤ θab (z k ) + γ2−i k ∇θab (z k )T d k ; Let τk = 2−i k ; Let z k+1 = z k + τk d k and k := k + 1; go to step 2.
We propose a different technique for obtaining a sufficient descent direction. If an inexact solution to the complementarity subproblem does not provide us with a sufficient descent direction, instead of immediately switching to a gradient direction, we continue with the solution of the subproblem. The structural characteristics of our subproblem solver provide us with such an option. If several attempts to solve the subproblem to increasing levels of exactness do not provide a sufficient direction of descent, we then switch to a gradient direction. The challenge in solving the subproblem arises from the massive sizes that may be a consequence of probability measure which while discrete could be specified over a large sample-space. To account for such a possibility, we employ a parallelizable splitting-based approach for solving the subproblem through
86
a solution of smaller scenario-based complementarity problems. Additionally, such a technique provides us with the ability to continue with the scheme if one does not obtain the required sufficient descent direction.
3.2 A Matrix-Splitting Subproblem Solver Matrix splitting techniques have been employed for solving large systems of equations as well as complementarity problems [19, 27]. In this subsection, we discuss the types of splitting schemes that may be employed for solving our subproblem of interest. Our subproblem can be cast as the following linear complementarity problem:
SubCP
0≤
à ! x
λ
⊥
à Hk
−∇c T
!Ã ! x
∇c
0
λ
+ Fk ≥ 0.
In constructing a splitting scheme for such a problem, we rely on showing the following matrix property for the resulting LCP [19]. Lemma 3.2.1 Suppose the Nash-Cournot problem satisfies assumptions 2.3.3 and 2.3.2. Then the matrix Hk as defined in equation 3.1 is a P-matrix. Proof : We note that Hk can be written as Hk = Gk + Rk , ∇211 f 1 . . . ∇21N f N . . Gk = . ∇2N 1 f N ... P λ1i ∇21 c 1i .. Rk = − .
∇2 f N
P
λN i ∇2N c N i
.
We claim that Gk is a P-matrix from theorem 2.3. In addition, Rk is a positive definite matrix, due to assumption 2.3.3. So Hk is a P-matrix. Lemma 3.2.2 Suppose the Nash-Cournot problem satisfies assumptions 2.3.3
87
and 2.3.2. Then the matrix M=
à Hk
−∇c T
∇c
²I
!
is row-sufficient. Proof : It may be recalled that a matrix M is row sufficient if [z i (M T z)i ≤ 0, ∀i = 1, . . . , N ] =⇒ [z i (M T z)i = 0, ∀i = 1, . . . , N ]. Let z=
à ! x
λ
,
Suppose z i (M T z)i ≤ 0, ∀i = 1, . . . , N . It follows that z T (M T z) ≤ 0; implying that
0≥
à !T à x HTk
λ
∇c T
!Ã ! x
−∇c
λ
= x T (HTk x). But Hk is P-matrix implying that if (x T HTk x) ≤ 0, then x ≡ 0. It follows that rowsufficiency of the original matrix holds if λi (−∇c T x)i = 0,
fori = 1, . . . , N .
But x = 0 implying the row-sufficiency of the original matrix. The row-sufficiency of the linear complementarity subproblem allows us to employ a splitting method. Such a method is given by the following set of steps. Suppose LCP(M , q) is the problem of concern. Let S = {z|z ≥ 0, M z + q ≥ 0}. The convergence of such a scheme follows from theorem 5.5.7 [19] and requires that M be row-sufficient, B be symmetric positive definite and LCP(q,M) be feasible. Theorem 3.2.3 (Th. 5.5.7 [19]) Let M be a row sufficient matrix and B be symmetric positive definite. Suppose that the LCP (q,M) is feasible. Then, the sequence {z k } produced by algorithm 5 is uniquely defined; moreover, every accumulation point of {z k } solves the LCP (q,M). If, in addition, the level set {z ∈ S : f (z) ≤ f (z 0 )} 88
Algorithm 5: A splitting method for LCPs 0
1 2 3
Let (B,C ) be a splitting of the matrix (M + M T )/2 with B being symmetric positive definite. Initialize τ > 0; Let z 0 ∈ S be an arbitrary starting point; Let k = 0; Given z k ∈ S, solve the quadratic program given by (QP k )
min r T z + 12 z T B z z ∈ S,
4 5
where r = q +C z k . Let z k+1/2 be the optimal solution of QPk ; let d k = z k+1/2 − z k ; if (d k )T Md k ≤ 0 then αk = 1 else let αk be the smallest number such that f (z k + αk d k ) = min{ f (z k + αd k ) : z k + αd k ∈ S, α ≥ 0}.
6 7
8 9
end Set z k+1 = z k + αk d k ; if (z k+1 )T (M z k+1 + q) < τ then stop else k = k + 1; go to step 3; end
is bounded, then so is the sequence {z k }. Proof : The aforementioned result is stated without a proof in [19]. While the splitting algorithm is more general in its applicability in that it lays weak requirements on the matrix properties, it has a key shortcoming that it requires a linesearch for globalization. Under a somewhat stronger requirement on the problem framework, one can obtain a markedly superior splitting scheme. We explain this next.
89
3.2.1 Linear constraints and contraction mappings Consider a simplified game-theoretic problem where agent i solves a problem of the form maximize p(X )x i − f i (x i ) + IEω [p ω (Y ω )y iω − h iω (y iω )]
Ai (x −i , y −i )
x i ,y i
subject to
x +C iω y iω ≤ C ap iω Aω i i x i , y iω ≥ 0,
ω ∈ Ω,
where p ω (Y ω ), p(X ) satisfy assumption 2.3.2 for all ω ∈ Ω. The functions f i , h i satisfy assumption 2.3.3. Then, the resulting complementarity subproblem is given by 0≤
x
Q
0 1 1 π A ⊥ .. . 0
y1
λ1 .. . yK λK
πK A K
0
−π1 A 1
π1 D 1
−π1 C1
π 1 C1
T
...
0
...
0
0
... .. .
0
0
0
...
π|Ω| D |Ω|
0
0
...
π|Ω| C|Ω|
T
−πK A K
T
x
1 y 1 0 λ .. .. . . |Ω| |Ω|T |Ω| −π C y 0
0
λ|Ω|
c
π1 h 1 π1 b 1 + .. . π|Ω| h |Ω| π|Ω| b |Ω|
where y 1ω Aω 1 .. ω ω y = . , A = y nω
..
.
C1 ω . , C = .. ω An 0
0 .. . , . . . Cn
... .. .
f 100 (x 1 ) . . . 0 .. − ∇2 (p(X )x), Q = . 0 0 x 0 . . . f N00 N (x N ) h 1ω00 (y 1ω ) . . . 0 .. − ∇2 ω (p ω (Y ω )y ω ) Dω = . 0 0 y i ω00 ω 0 . . . h N (y N )
∀ω ∈ Ω.
We assume for simplicity that the matrices D ω are positive definite. Since the size of complementarity problem may be arbitrarily large, we suggest a parallelizable scheme based on a splitting method. Unlike the earlier approach, in
90
≥0
this method, each iterate of the splitting scheme provides a steadily improving approximation to the solution of the true problem. This is in stark contrast with the earlier framework that uses the iterates to compute directions along which one determines an appropriate steplength by using a linesearch. Suppose the complementarity problem is given by 0 ≤ z ⊥ M z + q ≥ 0. Then, we use a splitting M = B + C where B is positive definite and scenariodecomposable. The resulting algorithm is given by Algorithm 6: A splitting-decomposition method for LCPs 0 Let (B,C ) be a splitting of the matrix M with B being symmetric positive definite, and scenario decomposable; 1 Initialize τ > 0; 2 k = 1, z k = z 0 ; k T k 3 while (z ) (M z + q) < τ do 4 q k = C z k + q; 5 for j = 1: |Ω| do 6
j
j
z k+1 ∈ SOL(B j , q k ) 7
end k = k + 1; end
Q 0 0 0 D 1 −C1T 0 C1 ²I 0 0 0 B = 0 0 0 0 0 0 0 0 0
0
0
0
...
0
0
0
0
...
0
0
0
0
...
0
0
...
0 0
2T
2
−C
C2
²I
0 ... .. .
0
0
0
. . . D |Ω|
0
0
0
. . . C|Ω|
D
91
0
0 Ã j 0 j =⇒ B = D 0 Cj |Ω|T −C ²I
0
−C j ²I
T
!
.
It follows that 12 (B + B T ) has Cholesky factors given by RQ 0 0 0 R = 0 0 0
0
0
0
0
0
...
0
R1
0
0
0
0
...
0
0
²2 I
0
0
0
...
0
0
0
R2
0
...
0
0 ... .. .
0
0
0 0 0 T , RQ RQ = Q, (R j )T R j = D j . 0 0 1 ²2 I
1
1 ²2 I
0
0
0
0
0
0
0
. . . R |Ω|
0
0
0
0
0
...
0
We invoke the following theorem from [19] to claim the convergence of the scheme. Theorem 3.2.4 (Theorem 5.3.12 from [19]) Suppose B is a positive definite matrix. Let Bb be a positive definite matrix such that BbT Bb = 1 (B + B T ). Suppose also 2
that ρ(Bb,C ) = kBb−T C Bb−1 k < 1. Then, for an arbitrary q and any starting z ≥ 0, the uniquely defined sequence of iterates {z k } converges to a unique solution of LCP((M + ²I ), q). In the light of this theorem, it is enough to construct (B,C ) such that the spectral condition is satisfied. We have assumed that all of RQ , R 1 , . . . , R |Ω| are positive definite. In this case the suggested splitting of B suffices. ρ(R,C ) = kBb−T C R −1 k ≤ kR −1 k2 kC k. The constraints of each problem (Ai (x −i , u −i )) can be scaled to obtain Ce = αC such that kCek <
1 kR −1 k2
.
As a consequence the resulting solution will need to be scaled accordingly.
92
3.3 Convergence Properties of the Algorithm In this section, we discuss the convergence schemes associated with the upperlevel Newton method as well as the splitting method employed for solving the subproblem. Proof of convergence is only provided for the exact Newton scheme. In future work, we plan to extend these ideas to the inexact scheme in algorithm 4. The global convergence of the exact version of algorithm 4 follows from the following result from [27] which we restate here for convenience. Theorem 3.3.1 (Theorem 10.4.15 [27]) Let F be twice continuously differentiable and b and a be given scalars with b > a > 0. Let {z k } be an infinite sequence generated by the exact version of algorithm 4. Then the following hold: 1. Every accumulation point is a stationary point of the merit function θab . 2. Let z ∗ be an accumulation point of {z k } such that (10.3.3) from [27] holds, then z ∗ is a solution to CP(F). 3. If {z k } is an isolated limit point , then the whole sequence converges to that point. 4. Suppose that z ∗ is a limit point of {z k } and a stable solution of VI(K,F). If p > 2 and γ < 12 then 1
(a) Eventually d k is always given by z k+ 2 − z k ; (b) Eventually a unit step is accepted so that z k+1 = z k + d k (c) The local convergence rate of Q-quadratic We intend to leverage the ideas of Josephy-Newton methods [27] to prove the convergence of the inexact scheme. The local convergence of the inexact scheme in algorithm 4 follows from theorem 7.3.8 in [27].
3.4 Numerical Results In this section we test the performance of algorithm 4 on a set of test problems. The results presented here are for the exact version of the Newton scheme. The test problems we used were created from the game arising between bidders in a
93
two period electricity market as depicted in section 2.2.1. The complementarity problem solved is of the form
0 ≤ x ⊥ F(x) ≥ 0,
x=
à ! x
λ
where F(x) = ∇φ(x) + M x + q, and we assume the problem to satisfy z i ∈ K i = {z ≥ 0|z ≤ C ap i } for each i . The function φ is taken from the following class of functions φ(x) = 12 x T H x + c T x + d + (1 +
X X (x(i ) − a i ))1/p + τ sin(x i )
H Â 0, p = 2, 4, 6 . . . , τ ≥ 0. φ(x) is strictly convex function iff τ = 0 and nonconvex otherwise. a i are chosen at random in [0, 1]. The problems were coded in M ATLAB/TOMLAB. All LCPs were decomposed into smaller LCPs (LLCPs). The number of LLCPs solved (κ) was also varied. We used KNITRO as our LLCP solver.
3.4.1 Scalability We first study scalability of the algorithm with the size of the probability support (i.e. the size of z). The test problem set is the same as the one mentioned in section 3.4. The results have generated by running the exact version of algorithm 4. The test problems were generated using different values for |N| and |Ω| and random seeds for φ(·). The number of variables in the problem is given the formula |V| = |N| + 3 × |N| × |Ω|. As said in the previous section, we also altered the number of LLCPs solved (κ) according to |V|. The results are presented in table 3.4.1. The field total
time denotes the total program time in seconds. It includes time taken for the splitting method, various function evaluations and the linesearches. The time taken in seconds by the splitting method alone , summed over all iterations in presented in field splitting time. The parallel time is estimated as
parallel time =
splitting time
94
κ
.
# 1 2 3 4 5 6 7 8
|N | 5 5 2 5 5 3 5 8
|Ω| 10 20 100 50 55 100 100 90
κ 31 61 86 151 166 301 215 271
|V| 155 305 602 755 830 903 1505 2168
splitting time
parallel time
total time
final residue
0.84 1.53 3.12 103.78 25.09 47.15 36.80 89.98
0.027 0.025 0.036 0.687 0.151 0.157 0.171 0.332
2.41 8.65 52.74 345.74 214.67 293.63 757.18 1824.60
4.70E-05 8.20E-05 9.40E-05 7.90E-05 7.00E-06 1.80E-05 3.10E-05 1.10E-05
Table 3.1: Scalability of algorithm 4 We see a large disparity between the time taken by the splitting method and the total program time, particularly when |V| is large. We have found that much of the time outside the splitting method is spent in function evaluations, particularly ∇φ and ∇2 φ. We see this as a limitation of our test problems, not of the algorithm. The splitting method scales remarkably well with the problem size. From table 3.4.1 (st ≡ splitting time) ½
¾
st1 /st2 max ≈ 6, |V|1 /|V|2 while for modest problems, viz. |V| ≤ 600 this ratio is ≈ 2. The total time (tt) shows better scaling: ¾
½
tt1 /tt2 ≈ 3. max |V|1 /|V|2 Expectedly the scaling improves when one assumes parallelizability. The parallel time (pt)
pt1 /pt2 max ≈ 1.2 |V|1 /|V|2 ½
¾
We next compare the scalability of the stand alone matrix splitting-decomposition method to emphasize the utility of algorithms 5 and 6. The results in table 3.4.1 were obtained by solving LCP of varying sizes. All times are in seconds. The same results are graphically represented in figure 3.4.1. It was observed that the splitting-decomposition method outperforms KNITRO thoroughly. The splitting method also scales well with |V| while KNITRO shows poor scalability.
95
# 1 2 3 4 5 6 7 8
|V| 10 20 29 30 41 51 61 76
knitro-time
splitting time
splitting residue
0.19 1.92 5.58 6.16 15.10 28.39 48.06 91.09
0.52 1.04 1.48 1.58 2.22 2.81 3.58 4.64
2.15E-06 3.53E-06 2.49E-06 2.67E-06 3.23E-06 4.05E-06 4.86E-06 5.22E-05
Table 3.2: Comparison of splitting-decomposition method and KNITRO
Figure 3.1: Comparison of splitting-decomposition method and KNITRO
3.5 Summary and Conclusions This chapter together with chapter 2 has extensively addressed the Nash-Cournot equilibrium problem. Chapter 2 was devoted mainly to enhancing the understanding of the problem. We motivated our model using some examples. We reviewed some game-theoretic solution concepts for stochastic games and identified the our stochastic Nash-Cournot solution as an S-adapted open loop equilibrium. Under fairly general conditions we then proved the existence and uniqueness of such an equilibrium. Chapter 3 provided a scalable and convergent algorithm for the solution of the problem introduced in chapter 2. The outer layer of the algorithm was constructed as an inexact Newton scheme, which iteratively solved linearizations of the original problem. Globalization was achieved by conducting a linesearch on the D-gap merit function. Stochasticity brought with it a the problem of scalability. Since no existing solver could reliably solve problems with sizes of our interest, we developed a splitting-decomposition solver for solving the LCP. 96
The convergence of both algorithms: the inexact Newton scheme and the LCP solver were proved. We tested the Newton scheme in its exact form on a collection of test problems, and also tested the splitting-decomposition LCP solver independently. Both methods were observed to show very good scalability. We conjecture that the performance of the inexact version of the solver would in fact be better than the results presented here for the exact version. Coding the inexact version and testing it is work under progress. Notably all results obtained were assuming serial processing of the LLCPs. Parallelization would further improve the performance.
97
References
[1] T. A LPCAN AND T. B ASAR, Distributed algorithms for nash equilibria of flow control games, Annals of Dynamic Games 7 (2003) to appear., (2003). [2] T. A LPCAN , T. B ASAR , R. S RIKANT, AND E. A LTMAN, Cdma uplink power control as a noncooperative game, Wireless Networks, 8 (November 2002.), pp. 659– 669. [3] M. A NITESCU, On solving mathematical programs with complementarity constraints as nonlinear programs, SIAM J. Optim., 15(4) (2005), pp. 1203–1236. [4] K. T. AU , J. L. H IGLE , AND S. S EN, Inexact subgradient methods with applications in stochastic programming., Math. Program., 63 (1994), pp. 65– 82. [5] J.-P. AUBIN AND H. F RANKOWSKA, Set-valued Analysis, Springer, 1990. [6] T. B ASAR AND G. O LSDER, Dynamic Noncooperative Game Theory, Classics in Applied Mathematics, SIAM, Philadelphia, 1999. [7] E. M. L. B EALE, On minimizing a convex function subject to linear inequalities, Journal of the Royal Statistical Society, 17B (1955), pp. 173–184. [8]
, On minimizing a convex function subject to linear inequalities, J. Roy. Statist. Soc. Ser. B., 17 (1955), pp. 173–184; discussion, 194–203. (Symposium on linear programming.).
[9] J. F. B ENDERS, Partitioning procedures for solving mixed-variables programming problems, Numer. Math., 4 (1962), pp. 238–252. [10] A. B ERKELAAR , C. D ERT, B. O LDENKAMP, AND S. Z HANG, A primal-dual decomposition-based interior point approach to two-stage stochastic linear programming, Oper. Res., 50 (2002), pp. 904–915. [11] D. B ERTSEKAS, Nonlinear Programming: 2nd Edition, Athena Scientific, Belmont, MA., 1999. [12] J. B IRGE, Decomposition and partitioning methods for multi-stage stochastic linear programs, Operations Research, 33 (1985), pp. 989–1007.
98
[13] J. R. B IRGE AND F. L OUVEAUX, Introduction to Stochastic Programming: Springer Series in Operations Research, Springer, 1997. [14] J. B LOMVALL, A multistage stochastic programming algorithm suitable for parallel computing, Parallel Comput., 29 (2003), pp. 431–445. Parallel computing in numerical optimization (Naples, 2001). [15] J. B LOMVALL AND P. O. L INDBERG, A Riccati-based primal interior point solver for multistage stochastic programming, EJOR, 143 (2002), pp. 452– 461. [16] P. B OGGS AND J. T OLLE, Sequential quadratic programming, Acta Numerica, 4 (1995), pp. 1– 50. [17] P. B OGGS , J. T OLLE , AND P. WANG, On the local convergence of quasinewton methods for constrained optimization, SIAM Journal on Control and Optimization, 20 (1982), pp. 161–171. [18] X. C HEN AND R. S. W OMERSLEY, Random test problems and parallel methods for quadratic programs and quadratic stochastic programs, Optim. Methods Softw., 13 (2000), pp. 275–306. [19] R. W. C OTTLE , J.-S. PANG , AND R. E. S TONE, The Linear Complementarity Problem, Academic Press, Inc., Boston, MA, 1992. [20] A. C OURNOT, 1897, Macmillan, Researches into the mathematical principles of the theory of wealth 1838, English Edition, ed. N. Bacon. [21] G. B. D ANTZIG, Linear programming under uncertainty, Management Sci., 1 (1955), pp. 197–206. [22] G. B. D ANTZIG AND P. G LYNN, Parallel processors for planning under uncertainty, Ann. Oper. Res., 22 (1989), pp. 1–21. [23] G. B. D ANTZIG AND G. I NFANGER, Large scale stochastic linear programs: Importance sampling and benders decomposition, tech. report, 1991. [24] J. E. DENNIS J R . AND J. J. MORÉ, A characterization of superlinear convergence and its application to quasi-Newton methods, Math. Comp., 28 (1974), pp. 549–560. [25] E. D. D OLAN AND J. J. M ORÉ, Benchmarking optimization software with performance profiles, Math. Program., 91 (2002), pp. 201–213. [26] W. S. D ORN, Duality in quadratic programming, Quart. Appl. Math., 18 (1960/1961), pp. 155–162.
99
[27] F. FACCHINEI AND J.-S. PANG, Finite Dimensional Variational Inequalities and Complementarity Problems: Vols I and II, Springer-Verlag, NY, Inc., 2003. [28] R. F LETCHER AND S. L EYFFER, User manual for filterSQP, May 21 1998. [29] P. E. G ILL , W. M URRAY, AND M. A. S AUNDERS, SNOPT: An SQP algorithm for large-scale constrained optimization, Tech. Report NA–97–2, San Diego, CA, 1997. [30] P. E. G ILL , W. M URRAY, M. A. S AUNDERS , AND M. H. W RIGHT, User’s guide for NPSOL (version 4.0): A Fortran package for nonlinear programming, Technical Report SOL 86-2, Department of Operations Research, Stanford University, Stanford, CA, USA, jan 1986. [31] P. E. G ILL , W. M URRAY, AND M. H. W RIGHT, Practical Optimization, Academic Press, Boston, MA, USA, 1981. [32] M. J. G OLDSMITH, Sequential Quadratic Programming Methods based on Indefinite Hessian Approximations, PhD thesis, Stanford University, 1999. [33] A. G ROTHEY AND K. M C K INNON, A superlinearly convergent trust region bundle method, tech. report, University of Edinburgh, Edinburgh, UK, 1998. [34] S. P. H AN, Superlinearly convergent variable metric algorithms for general nonlinear programming problems, Math. Programming, 11 (1976), pp. 263–282. [35]
, A globally convergent method for nonlinear programming, J. Optim. Theory Applic., 22 (1977), pp. 297–309.
[36] P. T. H ARKER, A variational inequality approach for the determination of oligopolistic market equilibrium, Math. Programming, 30 (1984), pp. 105– 111. [37] A. H AURIE , G. Z ACCOUR , AND Y. S MEERS, Stochastic equilibrium programming for dynamic oligopolistic markets, J. Optim. Theory Appl., 66 (1990), pp. 243–253. [38] J. H IGLE AND S. S EN, Stochastic decomposition: An algorithm for two stage linear programs with recourse, Mathematics of Operations Research, 16 (1991), pp. 650–669. [39]
, Stochastic Decomposition: A Statistical Method for Large Scale Stochastic Linear Programming, Kluwer Academic Publishers, Boston, MA., 1996.
100
[40] J.-B. H IRIART-U RRUTY AND C. L EMARÉCHAL, Convex Analysis and Minimization Algorithms, vol. 1 and 2, Springer, Berlin, 1993. [41] B. H OBBS, Linear complementarity models of nash-cournot competition in bilateral and poolco power markets, IEEE Transactions on Power Systems, 16 (2001), pp. 194–202. [42] B. H OBBS AND J.-S. PANG, Spatial oligopolistic equilibria with arbitrage, shared resources, and price function conjectures, Mathematical Programming, Series B, 101 (2004), pp. 57 – 94. [43] B. F. H OBBS , C. B. M ETZLER , AND J.-S. PANG, Strategic gaming analysis for electric power systems: An MPEC approach, IEEE Transactions on Power Systems, 15 (2000), pp. 638–645. [44] W. H OCK AND K. S CHITTKOWSKI, Test Examples for Nonlinear Programming Codes, Springer-Verlag New York, Inc., Secaucus, NJ, USA, 1981. [45] G. I NFANGER, Monte Carlo (importance) sampling within a Benders decomposition algorithms for stochastic linear programs, Ann. Oper. Res., 39 (1992), pp. 41–67. [46]
, Planning Under Uncertainty, Boyd and Fraser Publishing Co., 1994.
[47] A. S. J.L INDEROTH AND S. W RIGHT, The empirical behavior of sampling methods for stochastic programming, Optimization Technical Report 0201, Computer Sciences Department, University of Wisconsin-Madison, (2002). [48] C. R. F. J R . AND R. E. QUANDT, On the existence of cournot equilibrium, International Economic Review, 4 (1963), pp. 92–96. [49] R. K AMAT AND S. O REN, Two-settlement systems for electricity markets under network uncertainty and market power, Journal of Regulatory Economics, 25 (2004), pp. 5–37. [50] K. C. K IWIEL, Proximity control in bundle methods for convex nondifferentiable minimization, Mathematical Programming, 46 (1990), pp. 105– 122. [51] C. D. KOLSTAD AND L. M ATHIESEN, Necessary and sufficient conditions for uniqueness of a cournot equilibrium, Review of Economic Studies, 54 (1987), pp. 681–90. [52] C. D. KOLSTAD AND L. M ATHIESEN, Computing cournot-nash equilibria, Oper. Res., 39 (1991), pp. 739–748. [53] I. V. KONNOV, Equilibrium Models and Variational Inequalities, Elsevier, 2007. 101
[54] A. K ULKARNI AND U. V. S HANBHAG, On recourse properties of stochastic quadratic programs, to be submitted. [55] S. L EYFFER AND T. M UNSON, Solving multi-leader-follower games, Preprint ANL/MCS-P1243-0405, Argonne National Laboratory, Mathematics and Computer Science Division, (April, 2005). [56] J. L INDEROTH AND S. W RIGHT, Decomposition algorithms for stochastic programming on a computational grid, Comput. Optim. Appl., 24 (2003), pp. 207–250. Stochastic programming. [57] X. L IU AND G. Z HAO, A decomposition method based on sqp for a class of multistage stochastic nonlinear programs, SIAM J. on Optimization, 14 (2003), pp. 200–222. [58] A. L UCIA, An explicit quasi-newton update for sparse optimization calculations, Mathematics of Computation, 40 (1983), pp. 317–322. [59] Z.-Q. L UO, J.-S. PANG , AND D. R ALPH, Mathematical Programs with Equilibrium Constraints, Cambridge University Press, Cambridge, 1996. [60] E. M ASKIN AND J. T IROLE, A theory of dynamic oligopoly, i: Overview and quantity competition with large fixed costs, Econometrica, 56 (1988), pp. 549–69. [61]
, A theory of dynamic oligopoly, ii: Price competition, kinked demand curves, and edgeworth cycles, Econometrica, 56 (1988), pp. 571–99.
[62] C. M ETZLER , B. H OBBS , AND J.-S. PANG, Nash-cournot equilibria in power markets on a linearized dc network with arbitrage: Formulations and properties, Networks and Spatial Theory, 3 (2003), pp. 123–150. [63] F. H. M URPHY, H. D. S HERALI , AND A. L. S OYSTER, A mathematical programming approach for determining oligopolistic market equilibrium, Math. Programming, 24 (1982), pp. 92–106. [64] W. M URRAY, M. A. S AUNDERS , AND U. V. S HANBHAG, A second-order method for mathematical programs with complementarity constraints, To be submitted, (2006). [65] A. N AGURNEY AND D. Z HANG, Massively parallel computation of dynamic traffic networks modeled as projected dynamical systems, in Network optimization (Gainesville, FL, 1996), vol. 450 of Lecture Notes in Econom. and Math. Systems, Springer, Berlin, 1997, pp. 374–396. [66] J. F. N ASH, Equilibrium points in n-person games, Proceedings of National Academy of Science, (1950).
102
[67] J. V. N EUMANN AND O. M ORGENSTERN, Theory of Games and Economic Behavior, Princeton University Press, 1944. [68] J. N OCEDAL AND S. J. W RIGHT, Numerical Optimization, Springer Series in Operations Research, Springer-Verlag, New York, 1999. [69] W. N OVSHEK, On the existence of cournot equilibrium, Review of Economic Studies, 52 (1985), pp. 85–98. [70] K. O KUGUCHI, Expectations and stability in oligopoly models, in Lecture Notes in Economics and Mathematical Systems 138, Springer-Verlag, Berlin, 1976. [71] J. P ENG AND M. F UKUSHIMA, A hybrid newton method for solving the variational inequality problem via the d-gap function, 1997. [72] M. J. D. P OWELL, Algorithms for nonlinear constraints that use Lagrangian functions, Math. Prog., 14 (1978). [73]
, The convergence of variable metric methods for nonlinearly constrained optimization problems, in Nonlinear Programming 3, O. L. Mangasarian, R. R. Meyer, and S. M. Robinson, eds., Academic Press, New York, NY, 1978, pp. 27–63.
[74]
, A fast algorithm for nonlinearly constrained optimization calculations, in Lecture notes in Mathematics, Numerical Analysis, Dundee 1977, G. A. Watson, ed., Springer-Verlag, 1978, pp. 144–157.
[75] D. R ALPH AND S. W RIGHT, Some properties of regularization and penalization schemes for mpecs, http://www.optimization-online.org, (2003). [76] R. R OCKAFELLAR AND R.-B. W ETS, Scenarios and policy aggregation in optimization under uncertainty, Math. Operations Res. 16 (1991) 1-29, 16 (1991), pp. 1–29. [77] R. T. R OCKAFELLAR AND R. J.-B. W ETS, Stochastic convex programming: Kuhn-Tucker conditions, J. Math. Econom., 2 (1975), pp. 349–370. [78] A. RUSZCZYNSKI, Decomposition methods, in Handbook in Operations Research and Management Science, vol. 10, Elsevier Science, Amsterdam, 2003, pp. 141–212. [79] A. RUSZCZYNSKI AND A. S HAPIRO, Introduction, in Handbook in Operations Research and Management Science, vol. 10, Elsevier Science, Amsterdam, 2003, pp. 1–64. [80] A. RUSZCZYNSKI AND A. S HAPIRO, Stochastic programming models, in Handbook in Operations Research and Management Science, vol. 10, Elsevier Science, Amsterdam, 2003, pp. 1–63. 103
[81] U. V. S HANBHAG, Decomposition and Sampling Methods for Stochastic Equilibrium Problems, PhD thesis, Department of Management Science and Engineering (Operations Research), Stanford University, 2006. [82] U. V. S HANBHAG , G. I NFANGER , AND P. W. G LYNN, An iterative sampling method for stochastic quadratic programs, To be submitted, (2008). [83] D. F. S HANNO, On variable metric methods for sparse Hessians, Math. Comp., 34 (1980), pp. 499–514. [84] A. S HAPIRO, Monte carlo sampling methods, in Handbook in Operations Research and Management Science, vol. 10, Elsevier Science, Amsterdam, 2003, pp. 353–426. [85] H. D. S HERALI , A. L. S OYSTER , AND F. H. M URPHY, Stackelberg-NashCournot equilibria: characterizations and computations, Oper. Res., 31 (1983), pp. 253–276. [86] M. S LADE, What does an oligopoly maximize?, G.R.E.Q.A.M. 89a14, Universite Aix-Marseille III, 1989. [87] M. S PENCE, The implicit maximization of a function of monopolistically competitive markets. Discussion Paper No. 461, Harvard Institute of Economic Research, Cambridge, MA., 1976. [88] H. V. S TACKELBERG, The Theory of Market Economy, Oxford University Press, London, 1952. [89] J. S TOER AND R. A. TAPIA, On the characterization of q-superlinear convergence of quasi-Newton methods for constrained optimization, Mathematics of Computation, 49 (1987), pp. 581–584. [90] F. S ZIDAROVSZKY AND S. YAKOWITZ, A new proof of the existence and uniqueness of the cournot equilibrium, International Economic Review, 18 (1977), pp. 787–89. [91] P. L. T OINT, On sparse and symmetric matrix updating subject to a linear equation, Math. of Computation, 31 (1977), pp. 954–961. [92] R. M. VAN S LYKE AND R. W ETS, L-shaped linear programs with applications to optimal control and stochastic programming, SIAM J. Appl. Math., 17 (1969), pp. 638–663. [93] D. W. WALKUP AND R. J.-B. W ETS, Continuity of some convex-conevalued mappings, Proc. Amer. Math. Soc., 18 (1967), pp. 229–235. [94] D. W. WALKUP AND R. J.-B. W ETS, Stochastic programs with recourse, SIAM Journal of Applied Mathematics, 15 (1967), pp. 1299–1314. 104
[95] G. Y. W EINTRAUB , C. L. B ENKARD, AND B. V. R OY, Oblivious equilibrium: A mean field approximation for large-scale dynamic games., in NIPS, 2005. [96] R. B. W ILSON, A Simplicial Algorithm for Concave Programming, PhD thesis, Cambridge, 1963. [97] J. YAO, S. O REN , AND I. A DLER, Computing Cournot equilibria in two settlement electricity markets with transmission constraints., in Proceeding of the 38th Hawaii International Conference on Systems Sciences (HICSS 38). Big Island, Hawaii, 2004. [98] G. Z AKERI , A. B. P HILPOTT, AND D. M. RYAN, Inexact cuts in Benders decomposition, SIAM J. Optim., 10 (2000), pp. 643–657. [99] D. Z HANG AND U. V. S HANBHAG, Second-order barrier lagrangian schemes for stochastic nonlinear programming, working paper. [100] G. Z HAO, A log-barrier method with benders decomposition for solving two-stage stochastic linear programs, Mathematical Programming, 90 (2001), pp. 507–536. [101]
, A lagrangian dual method with self-concordant barriers for multistage stochastic convex programming, Math. Program., 102 (2005), pp. 1– 24.
[102] W. T. Z IEMBA, Computational algorithms for convex stochastic programs with simple recourse, Operations Research, 18 (1970), pp. 414–431.
105
Author’s Biography
Ankur A. Kulkarni was born in Mumbai, India on August 15th 1983. He is the elder of the two sons of Dr. Achyut S. Kulkarni and Mrs. Archana A. Kulkarni. He graduated from the Indian Institute of Technology, Bombay in 2006 with a B.Tech. in Aerospace Engineering. He has since been a graduate student at the University of Illinois at Urbana-Champaign with the Industrial and Enterprise Systems Engineering Department. He will be continuing his education at the same department by pursuing a Ph.D. in Industrial Engineering.
106
