Risk-neutral PDE-constrained generalized Nash equilibrium problems

A class of risk-neutral generalized Nash equilibrium problems is introduced in which the feasible strategy set of each player is subject to a common linear elliptic partial differential equation with random inputs. In addition, each player’s actions are taken from a bounded, closed, and convex set on the individual strategies and a bound constraint on the common state variable. Existence of Nash equilibria and first-order optimality conditions are derived by exploiting higher integrability and regularity of the random field state variables and a specially tailored constraint qualification for GNEPs with the assumed structure. A relaxation scheme based on the Moreau-Yosida approximation of the bound constraint is proposed, which ultimately leads to numerical algorithms for the individual player problems as well as the GNEP as a whole. The relaxation scheme is related to probability constraints and the viability of the proposed numerical algorithms are demonstrated via several examples.


Introduction
Whether it be a consequence of noisy measurements, estimated parameter values, or model ambiguity, uncertainty is present in just about every mathematical model of real-world phenomena. Whenever the uncertainty is untreatable by deterministic quantities, it is best to assimilate it into our mathematical models via random variables, vectors, or elements. This allows us to find more robust solutions in the face of future uncertainty and guard against outlier events. Since many models in engineering and the natural sciences are defined by partial differential equations (PDEs), the inclusion of random inputs leads us to consider parametric or random PDEs as part of our optimization problems, cf. [9,18,49,60,62].
PDE-constrained optimization under uncertainty is a challenging area of mathematical optimization with many relevant applications in the engineering sciences. It is a growing field with many recent of contributions in theory and algorithms, see e.g. [3, 11-14, 20, 22, 39-41, 43-45, 57, 58, 63]. However, many practical problems require the simultaneous minimization of multiple objectives. By pitting these objectives against each other, i.e., treating the problem as a noncooperative game with each objective and separate control representing a player and its individual strategy, we naturally come to study PDE-constrained Nash equilibrium problems under uncertainty. In the deterministic setting, we mention here the pioneering works [16,26,[51][52][53][54]61]. It is important to note, however, that the models in these papers do not consider bound constraints, in particular there are no state constraints. This is an important distinction, as it makes the difference between modeling the game via a coupled PDEs (no bound constraints) versus a variational inequality (no state constraints) versus a quasivariational inequality (with state constraints).
As with their deterministic counterparts, it is often necessary to look for a control that forces the state (solution of the PDE) to satisfy certain bound constraints, e.g., below a maximum temperature threshold or above a physical obstacle. When uncertain inputs are involved, this problem of state constraints becomes much more challenging. This is due in part to a lack of smoothness with respect to the random parameters and missing compactness properties, which we would expect in a deterministic setting. Moreover, though an adjoint equation solely for the state equation can be introduced, an adjoint equation that combines the state equation as well as a multiplier for the state constraint cannot be derived without assuming additional properties. The inclusion of state constraints leads in fact to generalized Nash equilibrium problems in Banach spaces. For recent work in the deterministic setting, we refer the reader to [33,34,37,38] and the references therein.
Summarizing the discussions above, we thus consider a class of risk-neutral PDEconstrained GNEPs under uncertainty subject to state constraints. In an abstract sense, this amounts to considering an N -player GNEP in which the ith player's problem takes the form Here, S(z) is the z-dependent random field solution of a linear elliptic PDE with uncertain inputs, Z i ad and K are closed convex sets and J i is an appropriate convex disutility function for player i. We will make the appropriate data assumptions below. The term "risk-neutral" arises due to the fact that only the expected disutility is considered. Letting z be a Nash equilibrium for this problem, player i would expect z i to be the best response to z −i on average, i.e., if the game were played repeatedly. Since the literature is rather scarce on the treatment of state constraints in PDE-constrained optimization under uncertainty, see e.g., [19,23] and the recent preprint [21], we pay special attention to the case where N = 1, as well. We comment further on the studies [19,23] below, which make use of probability constraints. In contrast, the abstract results in [22] can be used for state constraints as considered in this paper. However, these results require a different kind of constraint qualification that may be difficult to verify in general.
The contributions of our paper are as follows: 1. We exploit existing results on elliptic regularity theory to prove higher integrability and regularity of the random field solutions S(z). 2. Under appropriate constraint qualifications, we prove existence of solutions/ equilibria and derive optimality conditions for the optimization problem and GNEP. 3. We extend the well-known Moreau-Yosida approach for state constraints to the stochastic case and rigorously prove that the approximations converge to the original GNEP. 4. The link between the Moreau-Yosida regularization technique and probability constraints is established using concentration inequalities. 5. We propose and demonstrate the viability of numerical algorithms for the optimization problem and GNEP.
The first contribution is crucial, as we need at least essential boundedness of the random field solutions in order to use techniques of convex optimization in Banach spaces to develop the optimality theory. In (2), we require a Slater-type condition for the optimization problem and the strict uniform feasible response (SUFR) condition introduced in [33] for the GNEP. The SUFR condition imposes a kind of hidden symmetry on the GNEP model. Although Moreau-Yosida regularization has been used successfully in deterministic settings, the stochastic setting poses additional pitfalls. Nevertheless, passing to the limit in the relaxation parameter is crucial for the justification of the numerical methods in the fully continuous setting. The link to probability constraints in (4) is interesting in its own right, since the approximating problems are much easier to solve than a similar problem with probability conditions. In addition, we obtain a kind of probabilistic rate of convergence for the Moreau-Yosida relaxations, which is reflected in the properties of the out-of-sample controlled states in (5); even after solving with relatively small increasing batches and modest values of the relaxation parameter. The encouraging results in our numerical study (5) motivate a number of future research directions. The rest of the paper is structured as follows. In Sect. 2, we pose a number of basic assumptions along with an analysis of the forward problem. In addition, the optimization problems and GNEP are introduced. Following this, we derive existence and optimality conditions in Sect. 3; using the underlying structure and basic constraint qualifications. Due to the low multiplier regularity in the optimality conditions and a lack of adjoint equation in the sense that the righthand side is the sum of the derivative of the objective with respect to the state and the Lagrange multiplier for the state constraint, we propose a Moreau-Yosida technique in Sect. 4. This allows us to formulate function-space-based numerical algorithms for both the optimization problems and the GNEP in Sect. 5. The potential of the algorithms is demonstrated via several numerical examples. In particular, we provide a brief, post-optimal analysis using the performance of the computed controls to derive a statistic on the violation of the state constraint.

Notation, standing assumptions, and preliminary results
We start by defining the necessary function spaces. We assume that the physical domain D ⊂ R d with d = 1, 2, or 3 is an open bounded set such that D is either a convex polyhedron or the boundary of D, denoted by ∂ D, is of class C 1,1 .
The triple (Ω, F, P) denotes a complete probability space, where Ω is the sample space of possible outcomes, F the Borel σ -alegra of Ω for a fixed topology on Ω and P is a probability measure.
Given a real-valued Banach space (V , · V ), Borel measure μ, and p ∈ [1, ∞] we denote the usual Lebesgue-Bochner space L p μ (Ω; V ) of all strongly F-measurable V -valued functions by When V = R, we set L p μ (Ω; R) = L p μ (Ω) the usual Lebesgue space with underlying measure μ. When the Lebesgue measure μ = L is considered, we omit the subscript L and simply write L p (Ω). We denote by F L the σ -algebra of Lebesgue measurable sets. We recall here that for 1 ≤ p, q ≤ ∞ such that 1 /p + 1 /q = 1, it is known that For further information see [28,Chapter III]. We show in the sections below that the multipliers for the stochastic state constraints are of very low regularity, i.e., bounded additive measures. We will need the space ba, which we recall here for ease of reference, cf. [27,20.27 Definition] or [17].
The key result for our analysis related to this space is the existence of an isometric isomorphism between (L ∞ π (Ξ )) * and ba(Ξ , B, π), cf. [17,Thm. IV.8.16], where we use Finally, we fix several notational conventions. For a (real) Banach space V we denote the expectation of a random element X : Ω → V by For some nonempty subset C ⊂ V , I C : V → R ∪ {∞} represents the standard indicator function, which satisfies I C (x) = 0 if x ∈ C and +∞ otherwise. For an arbitrary convex set K , we define the standard convex normal cone by The (set-theoretic) characteristic function associated with some subset A is denoted by χ or χ A , where χ A (x) = 1 if x ∈ A and 0 otherwise. Strong convergence of a sequence is denoted by →, weak-convergence by , and weak-*-convergence by * . The closed ε-ball with center x in some normed space is denoted B ε (x). The superscript * is used to denote the adjoint operator or dual space. As usual C D means that C is bounded by D up to an independent constant. For two Banach spaces V and W , the set of all bounded linear operators from V to W will be denoted by L(V , W ). We use the typical convention from game theory for a vector u with N components for emphasizing the ith component by writing

PDE-constrained equilibrium problems as strategic games
As mentioned above, our results apply to both PDE-constrained optimization problems under uncertainty as well as stochastic equilibrium problems with PDE-constraints. Whereas the solution concept for PDE-constrained optimization is obvious, there are several possibilities for equilibrium problems from the perspective of game theory. The notation in this brief section is chosen to reflect the references to the game theory literature.
We recall that a strategic game comprises a set of N players or agents, their sets of actions A i , and a unique preference relation for each player over all possible profiles of actions a ∈ A := × N i=1 A i . In many cases, the preference relation can be described by the values of utility functions u i : A → R and the preferred solution concept for noncooperative behavior is often taken to be a Nash equilibrium; cf. [47]. The latter states thatā ∈ A is a (pure strategy) Nash equilibrium provided for all i = 1, . . . , N we have see, e.g., [48] for more details. We will refer to games in which the solution concept is a Nash equilibrium as Nash Equilibrium Problems or NEPs. We will take an analogous perspective for our PDE-constrained equilibrium problems. However, due to the presence of state constraints, the sets of actions are set-valued mappings A i (a −i ) that also depend on a −i for each i. This leads to a natural extension, first introduced by Debreu [15], see also [6]:ā ∈ A is a (generalized) Nash equilibrium provided for all i = 1, . . . , N we haveā i ∈ A i (ā −i ) and These games are significantly more difficult from both a theoretical as well as numerical perspective due to the embedded fixed point relation. We refer to games of this type as Generalized Nash Equilibrium Problems or GNEPs.

Linear elliptic random PDEs
Returning now to the context of PDE-constrained optimization, we introduce a class of linear elliptic random PDEs as our state system. Let Given z ∈ L 2 (D), we consider the following problem: Find u ∈ U such that for every ϕ ∈ H 1 0 (D) and consequently for every ϕ ∈ H 1 0 (D) . The reverse direction (from P-pointwise weak solutions to a solution of (2.3)) can be easily adapted from the nonlinear setting in [42]. The key components of the argument are: Prove the existence of a solution for P-a.e. ω, demonstrate measurability in ω using Fillipov's theorem for measurable selections, and obtain integrability using standard a priori estimates for elliptic PDEs. It is sometimes more convenient to work with one form versus the other as we will see below. For z = 0, we denote the solution of (2.3) by u f and for f ≡ 0 we set u = S(z). Hence, any solution u of (2.3) can be written We will demonstrate below that S(z) is a bounded linear operator in z between appropriate function spaces. In order to ensure well-defined solutions and derive higher regularity results, we make the following additional assumptions on the problem data.
Assumption 1 In addition to the standing assumptions on D, ∂ D, and (Ω, F, P), the following sets of assumptions will be necessary below.
The fixed bulk term f satisfies (ii) (Higher Regularity) In addition to (i), A ∈ L ∞ P (Ω; C 0,1 (D)). (iii) (Control Mapping) The control mapping B : Ω → L(L 2 (D) N , L 2 (D)) is measurable and essentially bounded, i.e. B ∈ L ∞ P (Ω, L(L 2 (D) N , L 2 (D)). Moreover, as a mapping from Ω to L(L 2 (D), H −1 (D)), B is completely continuous in the sense that for P-a.e. ω ∈ Ω we have Some remarks are in order. Assumption 1.(i) can be slightly weakened to allow for unbounded coefficients and still obtain the existence of solutions, cf. e.g., [24]. It is also possible to choose f and/or B(ω)z that is unbounded in ω. However, weakening these assumptions would mean that the solutions u to (2.3) are also not bounded. The latter property is essential for our treatment of state constraints. The Lipschitz continuity of A(ω, ·) :D → R in Assumption 1.(ii) will be used to ensure boundedness of u in x. This along with the regularity assumption on the boundary ∂ D can be slightly weakened to the extent that we can guarantee u ∈ L ∞ P×L (Ω × D), e.g., we could relax Lipschitz to Hölder and work with u(·, ω) in W 1, p (D) with p > d. The properties in Assumption 1.(iii) are the weakest possible for our analysis. Using Assumption 1, we gather several essential properties of the mapping z → u in the following result.

Proposition 1 Let Assumption 1 hold. For any z
) and the following a priori bound holds Here, C is independent of ω.
Proof Defining the bilinear form b : A(x, ·)∇u(x, ·) · ∇v(x, ·) dx and z-dependent linear form L(·; z) : U → R by we can view (2.3), as the variational problem: It readily follows from Assumption 1 that b is a U-coercive bilinear form. Then by the Lax-Milgram Lemma there exists a unique solution u ∈ U. In light of the equivalence to (2.5), we immediately deduce from the standard a priori bound: is P-essentially bounded. Due to the assumptions on A, C 1 does not depend on ω.
For the a priori bound (2.7), we need to consider two cases. We once again appeal to the equivalence between (2.3) and (2.5). If ∂ D is of type C 1,1 , then it follows from Assumption 1 along with Friedrichs' theorem, see e.g., [5,A12.2 Theorem], that for P-a.e. ω ∈ Ω we have . The same estimate also holds when ∂ D is nonsmooth, but D is a convex polyhedron, see Remark 1 below. The "constant" C(ω) is indeed a bounded and measurable function in ω. This follows from the fact that the term A(·, ω) C 0,1 (D) is measurable, uniformly bounded away from zero, and C(ω) is a sum of rational functions of A(·, ω) C 0,1 (D) , where it appears in a numerator and a denominator. Continuing, for P-a.e. ω ∈ Ω, we have Furthermore, we obtain Finally, due to the Gelfand triple +ess sup ω∈Ω f (ω, ·) L 2 (D) < ∞. Proposition 1 justifies the decomposition in (2.6). In particular, we see that S is a bounded linear operator and u f ∈ L ∞ P (Ω; H 1 0 (D) ∩ H 2 (D)). We deduce several additional properties in the following corollary. Proof Case (i) is a special case of [42,Prop 2.3]. In case (ii) linearity follows trivially from the definition of S(z) whereas boundedness is a consequence of (2.7) and Assumption 1.(iii).
We end this section by introducing a convenient P-pointwise notation that will aid in the derivation of optimality conditions below. We define to be the operators given by respectively. Note that A(ω) is a linear isomorphism due to the regularity results above. Given A, B we can understand S(z) + u f P-pointwise as whenever we need to work with higher regularity.

A class of risk-neutral PDE-constrained optimization problems
In this section, we introduce a class of optimization problems that will serve as a template for the individual player problems in the PDE-constrained GNEP.

Assumption 2
We assume that is a nonempty, closed, bounded, and convex set. (ii) (Objective) The cost parameter ν ≥ 0, u d ∈ L 2 (D), T ∈ L(L 2 (D)), and J : we define the state constraint by (iv) (Feasibility) There exists z ∈ Z ad such that (2.12) holds.
The boundedness in Assumption 2.(i) is only needed in the optimization setting if ν = 0. However, it is unclear how to extend the existence proof for the GNEP, as the latter follows from an application of the Kakutani-Fan-Glicksberg theorem, which includes a compactness condition. It is not necessary for our analysis to restrict ourselves to the tracking-type objective in Assumption 2.(ii). We could proceed in a more general manner as suggested in [44] under appropriate convexity, continuity, and growth conditions. This would require further technical assumptions that we believe would detract from the main purpose of the text. The nonemptiness of the feasible set in our setting is assumed in Assumption 2.(iv). Provided Z ad admits a z > 0 with sufficiently large L ∞ (D)-norm, then the existence of a feasible point can be guaranteed by the maximum principle in light of the regularity result in Proposition 1.
The inclusion of state constraints in PDE-constrained optimization in the form of (2.12) is new. An alternative way of interpreting (2.12) would be to consider either From the perspective of stochastic programming, this is rather restrictive and in general settings (beyond PDE-constrained optimization), may lead to empty feasible sets.
Typically one remedies this by selecting a minimum probability level p ∈ (0, 1) and considering instead: (2.13) Several recent studies have considered this perspective, see [19,23]. However, these approaches do not circumvent the fundamental difficulties encountered with state constraints in regards to multiplier regularity and mesh-independent numerical approaches. In addition, the functional is nontrivial to analyze and use in numerical algorithms. This usually requires P to admit a log-concave density and for S(z)(x, ω) to have a very specific structure with respect to ω. For more on probability constraints, we refer the reader to [50,59] and the related references therein.
We may now formulate the optimization problem min z∈Z ad (2.14)

A class of risk-neutral PDE-constrained GNEPs
We now introduce a noncooperative game with N players by using the results of the previous section. The individual ith player is assumed to solve the following optimization problem Here, the quantities Z i ad , T i , ν i , and u i d are defined analogously to those in the standard optimization setting, where we again require Assumptions 1 and 2 for each In what follows, we denote the collective admissible set of controls by Z ad = Z 1 ad × · · · × Z N ad . The main difference for the individual player problems lies in the definition of the control mapping B. For the sake of reference, we make the following assumption.

Assumption 3
The operator B has the additive representation where B i satisfies Assumption 1 for i = 1, . . . , N .
In light of the assumptions, we may also formulate the PDE-constrained GNEP in terms of the following reduced space problems.

Existence and optimality conditions
We first prove existence of optimal solutions of (2.14) and provide optimality conditions. Then, by extending the arguments used in [34], we prove the existence of generalized Nash equilibria for (2.15). Optimality conditions for a certain type of equilibria are also derived. We will use the concept of variational equilibria, which is strongly related to the notion of normalized equilibrium due to Rosen [56]; although Rosen's concept of normalized equilibrium was formulated using Lagrange multipliers. This is a specific class of Nash equilibria that can in many cases be computed numerically.

Risk-neutral PDE-constrained optimization problems
For the risk-neutral PDE-constrained optimization problems the existence and optimality conditions are formulated as follows.
Theorem 4 Let Assumptions 1 and 2 hold. Then (2.14) admits a solutionz. If ν > 0, thenz is unique. Moreover, if there exists a z 0 ∈ Z ad and a constant κ > 0 such that then there exists a measureμ ∈ ba(Ξ , B, π) such that is the continuous embedding.

(iii) (Subgradient Condition) The general inclusion holds
Here, the latter term must be understood for an arbitrary test function δz ∈ L 2 (D).

Remark 2
In this general setting, we cannot guarantee thatμ splits into a generalized product of measures that would allow us to write B * A − * ι * μ using an expectation. We explain this in more detail following the proof. However, the subgradient condition in (iii) can be brought into a form slightly more familiar to PDE-constrained optimization by introducing an adjoint variable λ that satisfies the pointwise adjoint equations Using this term we can "unfold" the general subgradient condition into the inclusion coupled to the adjoint Eq. (3.2), a state equation for Sz and complementarity conditions for the state constraint. Similarly, in some settings we may choose sufficiently regular test functions ϕ and introduce an additional adjoint variable λ to simplify the term Indeed,μ defines a bounded linear functional on L ∞ π (Ξ ). So if A can define a linear isomorphism between a subspace of L ∞ π (Ξ ) and, e.g., L 2 π (Ξ ), then λ would be a (very weak) solution. These two observations could potentially be used in a numerical setting, especially if P is discrete.
Proof To prove existence, we need to argue that the feasible set is weakly sequentially closed and F(z) := E[J (S(z) + u f , z)] is weakly sequentially lower semicontinuous on L 2 (D). Since the assumptions on J imply F is convex and the latter component of J is deterministic and continuous, we concentrate on the properties of S and their relation to the first argument of J .
By Assumption 2, (2.14) admits a feasible point and consequently a minimizing sequence {z k } ⊂ Z ad such that (2.12) holds. Since Z ad is bounded, closed, and convex, {z k } admits a weakly convergent subsequence {z k l }. For each l, we have Since S is completely continuous as a mapping into L 2 P (Ω; H 1 0 (D)), we have S(z k l ) → S(z) strongly. Moreover, the Sobolev embedding theorem (see e.g. [2, 4.12 Theorem]) and the fact that Continuing, the integrand J induces a superposition operator that is continuous from the product space , see e.g., [44,Ex. 3.2]. Then by combining the properties of S with this continuity result, we deduce the weak lower semicontinuity of F. It follows from the direct method thatz is an optimal solution, which is of course unique if ν > 0 as the objective would be strictly convex.
In order to derive first order optimality conditions for (2.14), we write min z∈Z ad and appeal to the general Lagrangian formalism in [10,Chap. 3]. Here, we set is the continuous embedding and K is the convex cone of all positive essentially bounded B-measurable functions. Note that we first use the continuous embedding of L ∞ P (Ω; ) and then the continuous embedding of L ∞ P (Ω; L ∞ (D)) into L ∞ π (Ξ ) to define ι. The latter two spaces are not equivalent.
Since K has a nonempty interior and G is clearly convex with respect to the partial order induced by (−K ), (3.1) is equivalent to the constraint qualification 0 ∈ int {G(Z ad ) − K } (and therefore Robinson's CQ), cf. [ Since K is a closed, convex cone,μ ∈ N K (G(z)) yields assertions (i) and (ii). To obtain the form in (iii), we first note that For the objective function F, we can exploit the equivalence with the pointwise adjoints and write This concludes the proof.
We caution the reader that the form of the duality pairing used for the μ-multiplier initially does not include the expectation with respect to P. However, if μ is σ -finite and σ -additive, then by the Radon-Nikodym theorem, there exists a density ρ μ such that In other words, we would have dμ = ρ μ d(L × P). Furthermore, The sign condition on μ carries over to ρ μ , in which case |ρ μ | = −ρ μ . This would indicate that ρ μ ∈ L 1 π (Ξ ). We could then write ). This would then allow us to incorporate the density into the adjoint equation, which is formulated in a very weak sense. This is essential, as otherwise the dual pairing with ρ μ and the test functions would not be defined.

Risk-neutral PDE-constrained GNEPs
In order to prove existence of at least one generalized Nash equilibrium and link the proof to a function-space-based numerical algorithm, we restrict ourselves to a variational reformulation as mentioned earlier. The variational reformulation is based on the so-called Nikaido-Isoda function Ψ : L 2 (D) N × L 2 (D) N → R. For our GNEP the Nikaido-Isoda function is given by We then introduce the potentially set-valued function R : Z ad ⇒ Z ad given by This mapping acts as a collective best-response function to a strategy vector z ∈ Z ad for all players simultaneously. Next, we define variational equilibria by their characterization as fixed points of the best-response function R . The nomenclature diverges somewhat from the literature, but it should be clear in context what is meant below.
Note that for jointly convex GNEPs, every variational equilibrium is also a Nash equilibrium [34,Theorem 3.2]. This characterization converts the proof of the existence of Nash Equilibria to a fixed point problem. The essential ingredient is the fixed point theorem of Kakutani-Fan-Glicksberg, see e.g. [4,Corollary 17.55].

Theorem 5 Let Assumptions 1 and 3 hold. The set of variational equilibria of the jointly convex GNEP (2.15) is weakly compact and nonempty.
Proof We proceed as in [34,Theorem 3.2], in order to apply the fixed point theorem of Kakutani-Fan-Glicksberg on R. By adapting the proof to the current setting, it follows from Theorem 4 that R has nonempty and convex images.
To ensure compactness, we recast the problem in the space X i , where X i is L 2 (D) endowed with the weak topology. Note that X i is a real locally convex topological space. The equivalence of weak and strong closure for convex sets in reflexive Banach spaces implies that Z i ad is closed in X i . Moreover, the weak compactness of closed and bounded convex subsets in reflexive Banach spaces implies that each set Z i ad is convex and compact in X i or equivalently sequentially compact (see [65, Satz VIII.6.1(Satz von Eberlein-Shmulyan)]). Consequently, if we take Z ad = Z 1 ad × · · · × Z N ad and X = X 1 ×· · ·× X N , then Z ad ⊂ X , where Z ad is also nonempty, convex and compact in X . Due to the latter property, the weak topology is metrizable on Z ad (see [65, Lemma VIII.6.2]).
In order to see the closedness of the graph of R, we introduce the set Now, we consider a closed subset C ⊂ X ad and a sequence {z n } n∈N ⊂ R −1 (C) with z n →z in X (i.e. z n z in L 2 (D) N ). For every z n we choose v n ∈ C ∩R(z n ). By a slight adaptation of the arguments in the proof of Theorem 4, we can show that X ad is sequentially compact. Hence, there exists a convergent subsequence v n k X →v with v ∈ C.
For some arbitrary w ∈ X ad it holds that By adapting the proof of Theorem 4, we can argue that This is a consequence of the properties of the expectation, the objectives J i and the solution operator S. In particular, it is essential that S is completely continuous into L 1 P (Ω; H 1 0 (D)). Using again the complete continuity, we have It follows thatv ∈ R(z), which proves the sequential closedness of the graph of R or equivalently the closedness in X ( [65, Theorem B.1.2]). We now apply Kakutani-Fan-Glicksberg's fixed point theorem. The set of Nash equilibria of the GNEP is nonempty and compact in X and thus, weakly compact in L 2 (D) N .
The optimality conditions for a generalized Nash equilibria reads as follow. We adapt the same notation as Theorem 4.
Proof Similiar to the proof of Theorem 4, we work with the general Lagrangian formalism. We first note thatz ∈ R(z). This is equivalent tō In order to derive first order optimality conditions for variational equilibria of (2.15), we recall some of the notation from the proof of Theorem 4 We again set and define the continuous embedding ι : In the notation of [10], we set which yields the parametric Lagrangian Since (3.4) is equivalent to the constraint qualification 0 ∈ int {G(Z ad )− K }, it follows from [10, Thm. 3.6] that Assertions (i) and (ii) are implied byμ ∈ N K (G(z i ,z −i )) since K is a closed, convex cone. To obtain the subgradient conditions in (iii), we first note that For the objective function, it holds that In order to see that the sum of the subdifferentials equals the product in (3.5), we refer to the proof of [34,Theorem 3.7]. Analogously to (3.3), we can write Moreover, [8, section 4.6] enables us to write the normal cones as For i = 1, . . . , N , we have the componentwise subgradient condition

A Moreau-Yosida regularization technique
The optimality conditions derived in Theorems 4 and 6 are not suitable for the development of algorithms. This is mainly due to the low regularity of the multiplierμ for the state constraint. To remedy this issue, we propose a Moreau-Yosida (MY) regularization technique, similar to the studies [1,30,31,34]. From the perspective of risk aversion, the MY-regularization can be seen as a measure of regret, e.g., as in [55], for the state constraint. We will also use several concentration inequalities below, which link MY-regularization to probability constraints. This further justifies the viability of the approach and provides a modeling solution for cases in which either the constraint qualification is hard to verify and/or it is not known if a feasible point for the original problem actually exists. To the best of our knowledge, this is the first time that such concentration inequalities have been used in the context of MY-regularization for infinite dimensional optimization under uncertainty.

Approximation of the risk-neutral PDE-constrained GNEPs
More specifically, the γ -dependent regularized problem of player i in the risk-neutral PDE-constrained GNEP (2.15) reads as where γ > 0. The usage of MY-regularization amounts to approximating the original GNEP by a more numerically tractable NEP. We will refer to this γ -dependent strategic game as NEP γ .

Existence and optimality conditions
The existence of a Nash equilibrium for every γ > 0 follows by using almost identical arguments to those in Theorem 5. Moreover, the first-order conditions have a similar, but numerically more workable form. We state the following theorem for ease of reference.
Theorem 7 Let Assumptions 1 and 3 hold. The set of variational equilibria of the jointly convex NEP γ (4.1) is weakly compact and nonempty. If z γ ∈ Z ad is a Nash equilibrium, then we have the following necessary and sufficient optimality conditions: Corollary 7 allows us to introduce the adjoint variablesλ γ i ∈ L 2 P (Ω; H 1 0 (D)) for i = 1, . . . , N and the associated adjoint equations:

Asymptotic considerations
We now investigate the behavior of NEP γ as γ → ∞. This is important for both theoretical as well as numerical considerations. We closely follow the approach in [33]. In order to ensure consistency of the relaxed problems, we will require the fulfillment of a constraint qualification as introduced in [33].

Definition 3
We say that (2.15) satisfies the strict uniform feasible response constraint qualification (SUFR), if there exists an ε > 0 such that for all i = 1, . . . , N and A few comments are in order. Traditional constraint qualifications such as the existence of a Slater point or in nonlinear programming the Mangasarian-Fromovitz constraint qualification (Robinson's CQ in infinite dimensions) were developed for optimization problems. They provide not only the existence of Lagrange multipliers, but also, they indicate a certain stability of the constraint set around the optimal solution. For example, the MFCQ gives us that the Lagrange multipliers associated with the point in question lie in a convex, compact polytope. In Theorem 4, it was enough to assume such a CQ without the need to adapt to the GNEP setting. However, for issues of approximation, we will see in the following that GNEPs require a much more robust CQ such as SUFR in order to exhibit the local stability needed to bound the dual variables; in this case the adjoint states and the constraint multipliers. From a game-theoretic perspective, we are requiring that each player has a feasible response to any strategy by its competitors such that the common state constraint is strictly uniformly fulfilled. Finally, as the current regularity assumptions on the random inputs only provide essential boundedness, we will need more regularity of the solutions.

Assumption 8 (Higher Parametric Regularity)
The set Ω is a compact Polish space. The solution mapping S(·) + u f is a continuous affine mapping from The need for Ω to be a compact Polish space will be evident in the following proof. Under weaker assumptions, we have already shown that the mapping S(·) + u f is a continuous affine mapping from H 2 (D)). The continuity assumption is actually weaker than it appears and can be guaranteed under mild assumptions (continuity in ω) on A(x, ω), B(ω) and f (x, ω), cf. the results in [35,Section 6]. The main idea is to reformulate the random PDE as a parametric fixed point equation and apply classic results on parametric dependence of solutions to fixed point equations. We now state the main result of this section. (4.2) as stated in Corollary 7. This sequence admits a limit point

Theorem 9 Suppose the GNEP (2.15) satisfies the Slater condition (3.1) and SUFR. If in addition Assumption 8 holds, then there exist sequences γ n → ∞ and
where, for all i = 1, . . . , N , we have Moreover, the limit point satisfies z * ∈ Z ad (4.5a) for an arbitrary test function ϕ ∈ L 2 (D). Finally, ρ * satisfies Note that (4.5c) and (4.5d) correspond to the subdifferential inclusion in Theorem 6. For readability, we split the proof over several partial results. Proof Fix a sequence γ n → ∞ for n → ∞. Since Z ad is weakly compact in L 2 (D) N and z γ n ∈ Z ad for all γ n , there exists a subsequence, denoted by γ k := γ n k and some element z * ∈ Z ad such that z γ k z * in L 2 (D) N . According to SUFR, there exists an ε > 0 and a sequence {v γ k } k→∞ ⊂ Z ad such that S(v Then for all γ k , the non-negativity of the MY-term gives us the lower bound: Furthermore, by definition of a Nash equilibrium we have the simple upper bound Using the fact that S is completely continuous into L 2 π (Ξ ) and each individual feasible set is bounded, we deduce the existence of a constant M independent of i, γ k such that Combining these observations yields Using the weak lower semicontinuity of the objective functions, it follows that the bound also holds for the limit We can conclude, that Thus, z * ∈ Z ad such that S(z * i , z * −i ) + u f ≥ ψ π-a.e., i.e. z * is a feasible strategy vector for the GNEP.
We note that for feasibility of z * , it is not necessary for ε to be positive in the SUFR condition. In what follows, we discuss the convergence of the stationary points individually. We start by showing that z * is also a generalized Nash equilibrium.

Lemma 2
Suppose the assumptions of Theorem 9 hold. Let {γ k } be the sequence of MY parameters from the proof of Lemma 1. Then there exists a subsequence {γ l } with γ l := γ k l → +∞ such that the weak limit point z * is a generalized Nash equilibrium.
s}. Due to the SUFR condition, X i is non-empty. Since for all γ k the associated z γ k is a Nash equilibrium, it holds that for all v i ∈ X i . For any v i ∈ X i , we want to construct a strongly convergent sequence Due to the SUFR condition, there exists an ε > 0 and for all k, a v k i ∈ Z i ad such that S(v k i , z γ k −i ) + u f ≥ ψ + ε, π -a.s. Clearly, {v k i } k→∞ is uniformly bounded in L 2 (D). Since every admissible set of each player is convex, we have that lies in Z i ad for all t ∈ (0, 1). Due to the linearity of the operator A and B, it holds that (4.8) We know, that for P-a.e. ω ∈ Ω the solution operator S(v i , ·)(ω) + u f (ω) maps continuously from L 2 (D) N −1 into H 1 0 (D) ∩ H 2 (D). Due to the Sobolev and Rellich-Kondrachov theorem, the solution of the state equation can be continuously and compactly embedded into the space of continuous functions overD P-a.s. Thus, S(v i , ·)(ω) + u f (ω) maps from L 2 (D) N −1 into C(D) for P-a.e. ω ∈ Ω. Combining this with the regularity assumption on the solution of the state equation, we have . Then by virtue of the nature of convergence in the C(Ω; C(D))-norm, we deduce the existence of a subsequence γ k l , denoted by γ l , such that then t l → 0 and t l ∈ (0, 1) for all l. Moreover, substituting (4.9) in (4.7) and due to (4.8), we have Passing to the limit as l → ∞ yields |t l | v l i L 2 (D) + v i L 2 (D) → 0 due to the boundedness of {v i (t l )} l→∞ and that {t l } l→∞ is a null sequence. Thus, we have constructed a sequence {v Finally, by substitution, we have For all i = 1, . . . , N , passing to the limit inferior yields the following inequality Here, we see that the uniformity in the SUFR condition is crucial to prove that z * is in fact a Nash equilibrium. In the following result, we obtain a stronger form of convergence to z * . This is necessary to derive the adjoint equation in the limit.

Lemma 3 Under the assumptions of Theorem 9, (4.4a) holds.
Proof First, we choose z * i ∈ X i in the construction of (4.7) with t = t l as in (4.9), then we have Then it holds that Passing to the limit superior yields lim sup Due to the complete continuity of S, we have Then (4.10) reads as This implies that lim sup Due to the weak convergence of {z Together with the weak convergence, the assertion follows.
We proceed with the sequence of the state variables.

Lemma 4 Under the assumptions of Theorem 9, (4.4b) and (4.5b) hold.
Proof This directly follows from the assumption, that S(·, ·) + u f : We note that the continuity in Ω is not really needed to prove a norm convergence result. Indeed since {z H 2 (D)). By Corollary 1 we even know that u γ ln → u * in L ∞ P (Ω; H 1 0 (D) ∩ H 2 (D)) holds. Next, we turn our attention to the sequence of the multipliers μ γ for the state constraint. We will observe that the Slater condition is enough to obtain a bound on

Lemma 5 Suppose the assumptions of Theorem 9 hold. In particular, (3.1) is fulfilled. Then we have (4.4c).
Proof We now prove the existence of a constant c 0 > 0 such that for any z ∈ B ε (0) ⊂ L ∞ π (Ξ ) and some fixed ε > 0. For the sake of readability, we set β : Unless otherwise noted, (·, ·) denotes the inner product on L 2 π (Ξ ) throughout the proof.
One readily shows that β is convex and continuously differentiable and therefore, μ γ = γβ (u γ ). Since β is convex, differentiable, and nonnegative, we obtain for any the equality By the assumption (3.4) there exists ε > 0 and z 0 ∈ Z ad such that for all v ∈ B ε (0) ⊂ L ∞ π (Ξ ): we have Since (Ω, F, P) is a complete probability space and the spatial domain D is bounded, the Lebesgue spaces are nested, and it holds that v ∈ L 2 π (Ξ ). Furthermore, Sz 0 +u f ∈ L 2 π (Ξ ). Fixing an arbitrary v ∈ B ε (0), we have Due to (4.13), we have The definition of the multiplier μ γ and the operator B yield Substituting the adjoint equation and applying the adjoint operator B * A − * yields Applying [28,Thm. 3.7.12] yields Here, the existence of c 0 is guaranteed, since the mappings and are continuously differentiable with uniformly bounded gradients on Z ad for all i = 1, . . . , N . This proves (4.12), since z was arbitrary. Using the fact that the L 1 -norm is positively homogeneous, subadditive and continuous, it follows from the Fenchel-Moreau theorem that the L 1 -norm is equivalent to the bidual norm It follows that the sequence {μ γ } γ →∞ is bounded in L 1 π (Ξ ). Therefore, by [17, Theorem IV.6.2] or [7, Corollary 2.4.3], we can extract a subsequence {μ γ l } l∈N which is weak * convergent to some regular countably additive Borel measure ρ ∈ M(Ξ). Next, we discuss the limit of the adjoint equation. We start by investigating the behavior of the expectation of the adjoint states. This leads to the derivation of a limiting adjoint state Λ * .
Then due to the Cauchy-Schwarz inequality and Hölder-inequality, respectively, we obtain Due to the continuous embedding of H 2 (D)∩ H 1 0 (D) into L 2 (D) and L ∞ (D), respectively, we have By the assumptions on the operators A and B, there exists C 2 ∈ L ∞ P (Ω) such that Now, combining the latter with (4.15) and (4.14), we obtain for all ϕ ∈ L 2 (D). Here, C 3 ∈ L ∞ P (Ω) and C 4 ∈ L 1 P (Ω). The existence of C 4 follows from the uniform bound on μ γ in the L 1 π (Ξ )-norm. Taking the expectation and applying Fubini's theorem yield for all ϕ ∈ L 2 (D). In other words, the sequence

Remark 3
The adjoint state plays an important role in numerical methods. In particular, P is often replaced by an empirical measure P N , which is associated with an i.i.d. random sample of size N . Therefore the quantity is of practical interest. By the (Kolmogorov) strong law of large numbers, we have with probability 1 as N → +∞ for any ϕ ∈ L 2 (D). For readability, set and recall that almost sure convergence implies convergence in probability. Then for fixed l ∈ N and any ε > 0, there exists N l,ε ∈ N such that On the other hand, the previous lemma gives us as l → +∞. It follows that for any ϕ ∈ L 2 (D) we have This means that the set of all events for which is contained in the set of all events for which Therefore, fix ϕ ∈ L 2 (D) and ε > 0, and choose l such that Then for all ε, there exists an l such that Thus, the diagonal sequence of sample averages of the adjoint variables weakly converges in probability to the limiting adjoint variable Λ * . For a fully discrete scheme using a finite element discretization of the underlying deterministic state spaces, in which error estimates for the deterministic adjoint variables were available, we could derive a similar statement. This is part of the justification for the update heuristic in our algorithm and, in general, for any related numerical algorithm in which the sample sizes gradually increase with the MY-parameters.
Next, we turn our attention on the adjoint equation in the limit.
Proof As in the previous proof, we start by constructing a specific test function. In this case, let w be the solution of the operator equation for all ϕ ∈ L 2 (D), then we know that w ∈ C(Ω; H 1 0 (D) ∩ H 2 (D)) holds. It follows that Taking the expectation on both sides yields We know that μ γ l * ρ * in M(Ξ). The right hand side reads as Passing to the limit l → ∞ yields for all ϕ ∈ L 2 (D).

Lemma 8 Under the assumptions of Theorem 9, (4.4e) and (4.5d) hold.
Proof Due (4.2), we can write Then the boundedness of the sequence {η Thus, there exists a η * i ∈ L 2 (D) and a subsequence {η γ ln i } n∈N such that the assertion holds.
Finally, we derive the complementarity system for the multiplier ρ * .
This completes the derivation of Theorem 9.

Probability constraints and Moreau-Yosida regularization
In this final theoretical section, we wish to draw the link between Moreau-Yosida regularization and probability constraints. We do so only for the the risk-neutral PDEconstrained optimization problem (2.14), as the treatment of the GNEP would require further technical assumptions and somewhat obfuscate our main point. The main tools are basic concentration inequalities from probability theory. We recall again the γdependent optimization problem: where γ > 0. We note that yet another way of formulating the original state constraint is Ideally, we would use the L ∞ (D)-norm as opposed to the L 2 (D)-norm, since the latter allows strong violation of the constraint on small subsets of positive measure for the weaker constraint for ε > 0, but arbitrarily small. However, in order to derive a result of the type in the following theorem with the L ∞ (D)-norm, we would need a careful analysis similar to [32]. This goes beyond the scope of the current paper.
Proposition 2 Let z γ be the unique minimizer of (4.16). Then for any ε > 0, we have and z is the unique minimizer of (2.14).
Proof Using Markov's inequality, we have We use z γ to obtain a simpler upper bound. By definition of z γ , it holds that for all v ∈ Z ad . In particular, we obtain the bound for all v ∈ Z ad such that S(v) + u f ≥ ψ for (L × P)-a.e. (x, ω) ∈ D × Ω. Using the minimizer z of (2.14) leads to From this we obtain E P (ψ − (S(z γ ) + u f )) + 2 L 2 (D) ≤ 2α γ . Then returning to Markov's inequality, we now have Finally, the complementary event is given by Remark 4 Using the analysis from the previous sections, we know that there exists a sequence γ n → +∞ such that the random variable converges strongly in L 1 (Ω, F, P) to Since z * is feasible, the state constraint holds and X * ≡ 0. Therefore, there exists a subsequence γ k := γ n k along which X k := X n k converges almost surely to 0; and consequently in distribution as well. For each k, we can set ε k = 1/ √ γ k and treat Y k := ε k as a degenerate random variable, which clearly converges in distribution to 0. It follows from Slutsky's theorem that X k + Y k converges in distribution to (ψ − (S(z * ) + u f )) + 2 L 2 (D) , i.e., 0 and since the Portmanteau lemma yields In this sense, Proposition 2 provides us with a probabilistic rate of convergence from Moreau-Yosida to feasibility for the original problem. We observe in the out-of-sample experiments in Sect. 5 almost exactly this behavior, i.e., for γ k = 1000, the percent of out-of-sample states is between one and three percent.

Numerical experiments
In this final section, we provide a numerical study to indicate how stochastic PDE-constrained optimization problems subject to pointwise state constraints and PDE-constrained GNEPs under uncertainty might best be solved. To the best of our knowledge, this is the first attempt to solve such problems numerically. As a result, the focus will be on the numerical solution of the individual optimization problems. For the GNEP, a Krasnoselskii-Mann-type alternating method is employed in which the dueling agents use the solver from Sect. 5.1.

Solving the individual problems
The basic idea behind this algorithm derives from the success of semismooth Newton methods for solving deterministic PDE-constrained optimization problems subject to state constraints using Moreau-Yosida regularization and path-following for the parameter updates; see e.g., [30,31]. Indeed, given γ > 0 and an iid sample of size M, we can replace the underlying probability distribution with the associated empirical probability measure P M and consider .
This is now a deterministic problem. In order to solve (5.1) with a semismooth Newton solver, we rewrite the first order optimality system as a single nonsmooth equation in z. The fixed random terms u f (ω m ) are defined analogously to u γ m . For readability, we denote the mapping z → B * λ γ as Λ(z) or Λ(z, ω) to indicate the dependence on ω. Moreover, we set In the current setting, where G is the Newton derivative of the projection operator. This allows us to apply a semismooth Newton method in L 2 (D) [29,64], which is known to be locally superlinearly convergent for each M and γ > 0.
However, since γ must be taken to +∞, such an algorithm would not be computationally efficient if M were chosen large for comparatively small γ . If M were to remain fixed, then we could use a strategy as in [1,30,31]. On the other hand, M should be ideally as large as possible or also treated as a parameter going to +∞. To remedy this issue, we set a maximum allowable sample size M max > 0 and penalty parameter γ max > 0 and, starting with M 0 ∈ N and γ 0 > 0, we add samples to M k every time γ k passes a certain threshold. For our numerical experiments, we consider a heuristic, which is motivated by the previous section; in particular the convergence statements in the fully continuous setting along with Remarks 3 and 4. A full convergence analysis linking sampling, approximation and smoothing error goes beyond the scope of this paper. The full algorithm is given in Algorithm 1. A few comments are in order.
The operator G is not explicitly given. Thus, it is necessary to use an iterative method to solve for the Newton steps dz k l , for which we use the tolerance tol newt ≥ 0. Since we are using a semismooth Newton iteration for pointwise bound constraints, the components of dz k l are fixed on the estimated active sets for each l and we only need to solve the linear systems on the potentially smaller inactive set. Here, it is important to note that each evaluation of G γ k M k (z k l )dz k l requires the solution of the forward equation and two adjoint equations for every sample m k = 1, . . . , M k . In our implementation, we employ a preconditioned conjugate gradient method. Therefore, the computational complexity of each Hessian-vector product involved must also be multiplied by M k and take into account the cost of applying the preconditioner. Similarly, the evaluation of the residual F γ M k (z k l ) requires a forward and adjoint solve for each sample. For our numerical examples, we use a direct solver for the linear elliptic PDEs.
Due to these facts, we suggest starting with a relatively small M 0 and increasing slowly with γ k . Moreover, we suggest a relatively large tol res 0 > 0 and ρ res close to 1. In step 13: of Algorithm 1, we simply set γ k+1 = φ(γ k ) = γ k + 1. More aggressive strategies may be possible, but empirical evidence suggests that this is not necessary and may even cause the Newton iteration to cycle. Finally, in step 15: of Algorithm 1, we link the increases of the sample sizes M k to γ k . For our implementation, we start with γ 0 and M 0 and increase M k by 10 every time γ k is divisible by 100. This is merely a heuristic and other strategies are possible.

Example: risk-neutral PDE-constrained optimization
In order to demonstrate the viability of the algorithm, we consider a model problem based on [43, Ex. 6.1, Ex. 6.2] and [35,Sec. 7.2]. Here, we set ν = 10 −3 , D = (0, 1), u(x) = sin(50.0 * x/π ), and consider the optimal control problem minimize z∈L 2 (D) where z ∈ Z ad := w ∈ L 2 (D) |−0.75 ≤ w(x) ≤ 0.75 a.e. x ∈ D and the solution of the random PDE u = u(z) ∈ L ∞ (Ω, F, P; H 1 (D)) solves the weak form of In addition, we impose the state constraint Furthermore, we suppose that with random variables ξ i : Ω → R, i = 1, 2, 3, 4, such that the supports ξ i , i = 1, 2, 3, 4, are [0, 1]. We assume here that each of these random variables is uniformly distributed. Following the usual change of variables, the forward problem (5.7) can be understood as with Ξ = [0, 1] 4 , endowed with the associated uniform density. We define ξ := (ξ 1 , . . . , ξ 4 ) ∈ Ξ . Since (5.8) is linear, we can use the superposition principle to lift the boundary conditions into the righthand side of (5.8). This allows us to transform the problem into the function space setting used throughout the paper.

Remark 5 (Feasibility)
In such settings as considered here, the ability to guarantee the nonnegativity of the state with at least one feasible control relies on three factors: 1. whether the maximum principle can be applied (a.s.), 2. whether the random inputs f , d 0 , and d 1 give rise to nonnegative solutions when solving the differential equations for each of these terms separately (a.s.), 3. the width of the bilateral bounds on z. If, for instance d 0 , d 1 are nonnegative, but f is bounded and negative on some portions of D with positive probability, then with sufficiently wide bounds on the control z, the righthand side can be made positive almost surely. This is admittedly not ideal. It is, however, much weaker than standard assumptions in two-stage stochastic programming such as complete or relatively complete recourse, which would require such a property for all feasible z.

Discretization and implementation
The pointwise forward problem and the control space are discretized using piecewise finite elements on a uniform mesh with parameter h = 1/(2 8 − 1). We use a standard Monte Carlo approximation for the random inputs ξ 1 , . . . , ξ 4 ∈ [0, 1]. We initialize the algorithm by choosing: γ 0 = 1, γ max = 10 4 , M 0 = 200, tol res 0 = 10 −2 , ρ res = 0.9997, tol newt = 10 −8 , z 0 ≡ 0. Once the penalty to sample threshold is reached, M k is increased by 10 samples. As mentioned above, the discrete PDEs are solved via a direct solver and the Newton steps are calculated using a preconditioned conjugate gradients method (for the linear equation on the inactive set). As a preconditioner we use the localized mass matrix for the inactive set. In the current implementation, we use the 2 -norm of the residual in the stopping criterion. Alternatively, one could use the proper discrete Riesz maps (i.e. the inverse mass matrix) to first obtain a representation of the discrete solution in the finite-dimensional subspace and then use the discrete L 2norm. This would be especially important in a nested grid or AFEM approach in future numerical studies. For the nonsmooth operator in the adjoint equations, we utilize a mass-lumping approach to obtain the discrete operators. Otherwise, the differential operators and identity operators give rise to the usual stiffness and mass matrices subject to the random inputs. We have observed that a heuristic damping strategy in which we take z k+1 = z k + t damp d k with t damp ∈ (0, 1), where z k is the current iteration and d k is the full Newton step, adds significant robustness to the method. For our experiments, we set t damp = 0.1. Though this has a clear effect on the local rate of convergence, the wildly varying ξ -dependent states appear to make this necessary.

Performance of the algorithm
The performance of the algorithm can be seen in Fig. 1, where we plot the total number of PCG iterations per γ -update (k) and the total number of Newton iterations needed to reach γ = 1000. The number of PCG iterations remains relatively stable (between 200 and 400), whereas the number of Newton iterations per γ -update appears to be trending downward. As mentioned above, robust convergence of the inner Newton iteration was ensured using a heuristic damping step with factor 0.1. Nevertheless, the number of iterations trends downwards as γ increases. We also note that both the CG algorithm and inner loop used inexact solves throughout. Since we employed a relatively rough initial stopping tolerance and a small batch of samples (despite increasing by 10 everytime (γ k mod 100) = 0) the algorithm consistently produces a solution z that performs exactly as expected in light of the model (risk neutral objective) and theory (especially Remark 4). This is qualitatively illustrated in Fig.  2, where we observed that only 0.6% of the out-of-sample states violated the bound constraint. Due to the presence of the random viscosity term in the forward problem, the L ∞ -norm of the sampled states can vary significantly. Finally, and perhaps due to the previous fact, we noticed that smaller batches sizes, e.g., on the order of 10, led to a failure of the Newton solver for γ near 1000.

A general algorithm
As mentioned earlier, we employ a fixed point strategy to solve a two-player, risknetural PDE-constrained GNEP. The fixed point iteration is derived from a standard Krasnoselskii-Mann iteration. We introduce the mappings T i (z j ) i = j, where and T 2 (z 1 ) is defined analogously. The fixed point iteration is based on the following outer iteration: 1. Given (z old 1 , z old 2 ) ∈ L 2 (D) × L 2 (D). 2. The first player determines z 1 = T 1 (z old 2 ) and reveals this to the second player. 3. The second player then determines z 2 = T 2 ( z 1 ) and reveals this to the first player. 4. The first player now determines w 1 = T 1 ( z 2 ). Obviously, (1)-(5) represents an ideal setting as the state constraint needs to be treated by a Moreau-Yosida approximation. In this context, we denote the γ -dependent mapping in steps (2)-(4) by T i γ . The full algorithm is depicted in Algorithm 2. Note that λ = 1 would correspond to a Gauss-Seidel iteration and λ > 1 to successive overrelation. As the evaluation of the T-mappings requires an iterative solver in practice, we suggest to initialize the nonlinear solvers by using z old 1 is used in (2), z old 2 in (3), z 1 in (4), and z 2 in (5).
Many of the inputs in Algorithm 2 are either self-explanatory or play the same role as in Algorithm 1. Here, we introduce tol km 0 > 0 and ρ km ∈ (0, 1], which allow us to successively reduce the tolerance used in the Krasnoselkskii-Mann iteration as γ k (and consequently M k ) increase. We suppress the fact that certain fixed data and parameter values need to be passed to the T γ -operators throughout the inner iterations.
It is again possible to adapt the tolerance used in the PCG solver for the Newton steps, but empirical evidence indicates that this value should be rather small (order at least 1e-6). Though the structure of Algorithm 2 is very similar to that of Algorithm 1, it is important to note that each evaluation of T i γ k is associated with a semismooth Newton solve for the current γ k and sample of size M k .

Remark 6
For each fixed γ k and M k , the algorithm is basically a Krasnoselskii-Mann iteration with inexact evaluations of the fixed point mapping. As such, convergence can be guaranteed if the latter can be shown to be nonexpansive. Such an analysis goes beyond the scope of the paper. Given the underlying individual problems are strongly convex, this property is most likely linked to the modulus of strong convexity of the individual cost functions.

Remark 7
Gauss-Seidel iterations (using λ = 1 above) have been considered in the context of GNEPs and NEPs in a number of texts, e.g., [33,34], for more than two players. A similar adaptation using the Krasnoselskii-Mann-type approach for more than two players is therefore conceivable.

Discretization and implementation
The discretization, sampling, γ -update strategy, and tolerance reduction for the Newton iterations are the same as in Sect. 5.2.1. We fixed λ = 0.5. Though further experiments demonstrated that successive over relaxation, i.e., λ > 1, does in fact work, the number of Krasnoselskii-Mann iterations remained roughly the same. The inner KM-iterations stopped once the discrete L 2 (D)-norm of z new 1 − z old 1 reached a tolerance of 1e-3.

Performance of the Algorithm
We have already investigated the performance of the sample-average based semismooth Newton solver in Sect. 5.2.2. As expected, the algorithm performs reliably in the GNEP setting, where it is called hundreds of times without failing to converge. The behavior of Algorithm 2 for the full control action is depicted in Figs. 3 and 4. For the second example, in which the controls are restricted to subsets of D, we point the reader to Figs. 5 and 6. In both cases, we observe non-trivial active sets for the equilibrium controls. In either case, the fixed point iteration requires a moderate number of iterations for the first few γ . This then rapidly tapers off as γ k and M k increased. The performance of the equilibrium controls is also demonstrated in Figs. 4 and 6. These plots correspond to an estimated violation (in the sense of the L ∞ -norm) of the state constraints of 1.5% and 0.6% respectively. This is well within the usual tolerance of 95% often used for probability constraints.

Conclusions and outlook
In this paper, we proved existence of solutions/equilibria and derive optimality conditions for both stochastic PDE-constrained optimization and equilibrium problems subject to state constraints. For our analysis, higher regularity of the random states was proven using a priori estimates for deterministic elliptic PDE. This allowed us  to make use of the existing optimality theory for convex optimization problems. In the case of GNEPs, a GNEP-specific constraint qualification was crucial for the development of a relaxation approach on which both the theory and our numerical methods could be built. We saw that this condition is fundamentally different than the classical constraint qualifications from nonlinear programming such as the Mangasarian-Fromowitz CQ, which was originally introduced in [46]. Nevertheless, the low regularity of the Lagrange multipliers still makes passing to the limit highly nontrivial.
After rigorously passing to the limit in the smoothing parameter, we provided further insight into the approximation technique using results on concentration inequalities and asymptotic statistics. Finally, we suggested two algorithms; the first for solving risk-neutral PDE-constrained optimization problems subject to state constraints and the second for the extension to GNEPs.
The algorithms performed well and the statistical properties of the solutions are comparable with what one would require of probability constraints; though our approach is much easier to treat theoretically and numerically. At least for a fixed sample, the optimization solver is known to converge locally superlinearly. A full convergence analysis linking sampling, adaptive finite elements, smoothing, and convergence of these algorithms (as least in a probabilistic sense) will be a future direction of research. The convergence of the GNEP solver is much more delicate and will require a fine analysis of the nonexpansivity of the underlying fixed point mapping. We postulate here that this is linked to the modulus of strong convexity of the underlying problems.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.