FOR

. We present a multigrid algorithm to solve eﬃciently the large saddle-point systems of equations that typically arise in PDE-constrained optimization under uncertainty. The algorithm is based on a collective smoother that at each iteration sweeps over the nodes of the computational mesh, and solves a reduced saddle-point system whose size depends on the number N of samples used to discretized the probability space. We show that this reduced system can be solved with optimal O ( N ) complexity. Wetestthemultigrid method on three problems: a linear-quadratic problem for which the multigrid method is used to solve directly the linear optimality system; a nonsmooth problem with box constraints and L 1 -norm penalization on the control, in which the multigrid scheme is used within a semismooth Newton iteration; a risk-adverse problem with the smoothed CVaR risk measure where the multigrid method is called within a preconditioned Newton iteration. In all cases, the multigrid algorithm exhibits very good performances and robustness with respect to all parameters of interest.


Introduction.
In this work, we present a multigrid method to solve the saddle point system where x = (y, u, p) = (y 1 , . . . , y N , u, p 1 , . . . , p N ) ⊤ , S has the block structure and all submatrices involved represent the discretization of some differential operators. More details on each block are provided in Section 2. Matrices such as (1.2) are often encountered while solving PDE-constrained optimization problems under uncertainty of the form s.t. y(ω) ∈ V satisfies e(y(ω), u, ω), v = 0 ∀v ∈ V, a.e. ω ∈ Ω, where u is the unknown deterministic distributed control, y(ω) is the state variable which satisfies a random PDE constraint expressed by e(·, ·, ω) for almost every realization ω of the randomness, Q is a real-valued quantity of interest (cost functional) and R is a risk measure. The vectors {y j } N j=1 and {p j } N j=1 are the discretizations of the state and adjoint variables y(ω) and p(ω) at the N samples in which the random PDE constraint is collocated. The vector u is the discretization of the deterministic control u. Problems of the form (1.3) are increasingly employed in applications. The PDE constraints typically represent some underlying physical model whose behaviour should be optimally controlled, and the randomness in the PDE allows one to take into account the intrisinc variability or lack of knowledge on some parameters entering the model. The introduction of a risk measure in (1.3) allows one to construct robust controls that take into account the distribution of the cost over all possible realizations of the random parameters. Therefore, the topic has received a lot of attention in the last years, see, e.g. [19,20,28,16,15,1,30,13,2].
However, few works have focused on efficient solvers for the optimality systems (1.1). A popular approach is to perform a Schur complement on u and solve the reduced system with a Krylov method (possibly with Conjugate Gradient), despite each iteration would then require the solution of 2N PDEs, with A j and A ⊤ j for j = 1, . . . , N [18]. For a full-space formulation, block diagonal preconditioners have been studied in [29], using both an algebraic approach based on Schur complement approximations, and an operator preconditioning framework.
In this manuscript, we present a multigrid method to solve general problems of the form (1.1) and show how this strategy can be used for the efficient solution of three different Optimal Control Problems Under Uncertainty (OCPUU). First, we consider a linear-quadratic OCPUU and use the multigrid algorithm directly to solve the linear optimality system. Second, we consider a nonsmooth OCPUU with box constraints and L 1 regularization on the control. To solve such problem, we use the collective multigrid method as an inner solver within an outer semismooth Newton iteration. Incidentally, we show that the theory developed for the deterministic OCPs with L 1 regularization can be naturally extended to the class of OCPUU considered here. Third, we study a risk-adverse OCPUU involving the smoothed Conditional Value at Risk (CVaR) and test the performance of the multigrid scheme in the context of a nonlinear preconditioned Newton method.
The multigrid algorithm is based on a collective smoother [5,6,38] that, at each iteration, loops over all nodes of the computational mesh (possibly in parallel), collects all the degrees of freedom related to a node, and updates them collectively by solving a reduced saddle-point problem. For classical (deterministic) PDE-constrained optimization problems, this reduced system has size 3 × 3 so that its solution is immediate [6]. In our context, the reduced problem has size (2N + 1) × (2N + 1), so that it can be large when dealing with a large number of samples. Fortunately, we show that it can be solved with optimal O(N ) complexity. Let us remark that collective multigrid strategies have been applied to OCPUU in [7,4] and in [34]. This manuscript differs from the mentioned works since, on the one hand, [7,4] considers a stochastic control u, therefore for (almost) every realization of the random parameters a different control u(ω) is computed through the solution of a standard deterministic OCP. On the other hand, [34] considers a Stochastic Galerkin discretization, and hence the correspoding optimality system has a structure which is very different from (1.2).
The multigrid algorithm assumes that all variables are discretized on the same finite element mesh. Further, despite we here consider only the case of a global distributed control, the smoothing procedure could be extended to a local distributed control on a subset D of the computational domain D, and to boundary control, see, e.g., [5].
The rest of the manuscript is organized as follows. In Section 2 we introduce the notation, a classical linear-quadratic OCPUU, and interpret (1.2) as the matrix associated to an optimality system of the discretized OCPUU. In Section 3, we present the collective multigrid algorithm, discuss implementation details and show its performance to solve the linear-quadratic OCPUU. In Section 4, we consider a nonsmooth OCPUU with box constraints and a L 1 regularization on the control. Section 5 deals with a risk-adverse OCPUU. For each of these cases, we first show how the multigrid approach can be integrated in the solution process, by detailing concrete algorithms, and then we present extensive numerical experiments to show the efficiency of the proposed framework. Finally, we draw our conclusions in Section 6.

2.
A linear-quadratic optimal control problem under uncertainty. Let D be a Lipschitz bounded domain in R d , V the standard Sobolev space H 1 0 (D), and (Ω, F , P) a complete probability space.
Given a function u ∈ L 2 (D), we consider the linear elliptic random PDE where (·, ·) denotes the standard L 2 (D) scalar product. To assure uniqueness and sufficient integrability of the solution of (2.1), we make the following additional assumption.
Assumption 2.1. There exist two random variables a min (ω) and a max (ω) such that . ω ∈ Ω, and further a min and a max are in L p (Ω) for some p ≥ 4.
Under Assumption 2.1, it is well-known (see, e.g., [25,10]) that (2.1) admits a solution in V for P-a.e. ω, and the solution y, interpreted as a function-valued random variable y : ω ∈ Ω → y(ω) ∈ V , lies in the Bochner space L p (Ω; V ) [12]. We often use the shorthand notation y ω = y(x, ω) when the dependence on x is not needed, or y ω (u) when we want to highlight the dependence on the control function u.
In this manuscript, we consider the minimization of functionals constrained by (2.1). Let us first focus on the linear-quadratic problem where y d ∈ L 2 (D) is a target state, f ∈ L 2 (D) (we omit the continuous embedding operators from L 2 (D) to L 2 (Ω; L 2 (D))), E : L 1 (Ω) → R is the expectation operator and ν > 0. Introducing the linear control-to-state map S : u ∈ L 2 (D) → y ω (u) ∈ L 2 (Ω; L 2 (D)), we reformulate (2.2) as, Existence and uniqueness of the minimizer of (2.3) follows directly from standard variational arguments [24,17,39,19]. Furthermore, due to Assumption 2.1, the optimal control u satisfies the variational equality where the adjoint operator S ⋆ : L 2 (Ω; L 2 (D)) → L 2 (D), satisfying (Su, z) L 2 (Ω;L 2 (D)) = (u, S ⋆ z) L 2 (D) , is characterized by S ⋆ z = E [p] where p = p ω (x) is the solution of the adjoint equation The optimality condition (2.4) can thus be formulated as the optimality system To solve numerically (2.2), we replace the exact expectation operator E of the objective functional by a quadrature formula E with N nodes Common quadrature formulae are Monte Carlo, Quasi-Monte Carlo and Gaussian formulae. The latter requires that the probability space can be parametrized by a sequence (finite or countable) of random variables {χ j } j , each with distribution µ j , and the existence of a complete basis of tensorized L 2 µj -orthonormal polynomials. Hence for the semi-discrete OCP, the P-a.e. PDE-constraint is naturally collocated onto the nodes of the quadrature formula.
Concerning the space domain, we consider a family of regular triangulations {T h } h>0 of D, and a Galerkin projection onto the finite element space V h of continuous piecewise polynomial functions of degree r over T h that vanish on ∂D, spanned by a Lagragian basis. N h is the dimension of V h , and {φ} N h i=1 is a nodal Lagrangian basis. We discretize the state, adjoint and control variables on the same finite element space.
Once fully discretized, (2.6) can be expressed as where A j are the stiffness matrices associated to the bilinear forms a ωj (·, ·), M s is the mass matrix, y d and f are the finite element discretizations of y d and f respectively, while y j and p j are the discretizations of y ωj and p ωj .
3. Collective multigrid scheme. In this section, we describe the multigrid algorithm to solve the full space optimality system (2.7). For the sake of generality, we consider the more general matrix (1.2), so that our discussion covers also the different saddle-point matrices obtained in Sections 4 and 5.
For each node of the triangulation, let us introduce the vectors y i and p i , which collect the degrees of freedom associated to the i-th node, the scalar u i = (u) i , and the restriction operators The prolongation operators are P i := R ⊤ i , while the reduced matrices S i := R i SP i ∈ R 2N +1×2N +1 represent a condensed saddle-point matrix on the i-th node, and satisfy Given an initial vector x 0 , a Jacobi-type collective smoothing iteration computes for n = 1, . . . , n 1 , where θ ∈ (0, 1] is a damping parameter. Gauss-Seidel variants can straightforwardly be defined. Next, we consider a sequence of meshes {T h ℓ } ℓmax ℓ=ℓmin , which we assume for simplicity to be nested, and restriction and prolongator operators R ℓ ℓ−1 , P ℓ ℓ−1 which map between grids T h ℓ−1 and T h ℓ . In the numerical experiments, the coarse matrices are defined recursively in a Galerkin fashion starting from the finest one, Nevertheless it is obviously possible to define S ℓ as the discretization of the continuous saddle-point system onto the mesh T h ℓ . With this notation, the V-cycle collective multigrid is described by Algorithm 3.1, which can be repeated until a certain stopping criterium is satisfied. We used the notation Collective Smoothing(·, ·, ·) to denote possible variants of (3.2) (e.g. Gauss-Seidel).
Notice that (3.2) requires to invert the matrices S i for each computational node. We now show that this can be done with optimal O(N ) complexity. Indeed, performing a Schur complement on u i , the system can be solved exclusively computing inverses of diagonal matrices and scalar products between vectors through (3.3) x n1 =Collective Smoothing(x 0 , S ℓ , n 1 ) (n 1 steps of collective smoothing) 5: (compute the residual) 6: 7: 8: x n2 =Collective Smoothing(x 0 , S ℓ , n 2 ) (n 2 steps of collective smoothing) 9: Notice that we should guaranteee that diag(a i ) admits an inverse and that ( . This has to be verified case by case, so we now focus on the specific matrix (2.7). On the one hand, the vectors a i are strictly positive componentwise, since (a Assuming that the V-cycle algorithm requires a constant number of iterations to reach a given tolerance as both N h , N and the number of levels increase (so that the size of the coarse problem becomes neglegible in the limit), Algorithm 3.1 exhibits an overall optimal linear complexity O(N h N ) since asymptotically the cost of each iteration is dominated by the N h relaxation steps, each of a cost proportional to N due to the above discussion. In the next numerical experiments sections, we show indeed that the number of iterations remain constant for several test cases.
The multigrid algorithm has a perfect robustness with respect to all parameters considered. In particular, the convergence does not deteriorate as ν → 0, which is a well-known troublesome limit for several preconditioners (see, e.g., [29,32,31,11] and references therein). Further, the third sub-table shows the robustness of the algorithm with respect to the number of levels as the fine grid is refined. The fourth sub-table instead shows that the number of iteration remains constant as the discretization of the probability space is refined. As a complementary experiment, we next consider L 2 = 0.025 and use a Monte Carlo quadrature formula, since Stochastic Collocation suffers from the curse of dimensionality. We also refine the spatial grid to properly capture the variations of the random diffusion coefficient. Table 2 confirms the robustness of the algorithm even for much larger numbers of samples.
4. An optimal control problem under uncertainty with box-constraints and L 1 penalization. In this section, we consider the nonsmooth OCPUU where y ω (u) ∈ V solves a ω (y ω (u), v) = (u + f, v), ∀v ∈ V, P-a-e. ω ∈ Ω,  (14) 21 (14) N h = 3969, ν = 10 −4 , N L = 3, σ 2 = 1.5, L 2 = 0.025. with a < 0 < b and ν, β > 0. Deterministic OCPs with a L 1 penalization lead to optimal controls which are sparse, i.e. they are nonzero only on certain regions of the domain D [37,8]. Sparse controls can be of great interest in applications, because it is often not desirable, or even impossible, to control the system over the whole domain D. For sparse OCPUU, we mention [23] where the authors considered both a simplified version of (4.1) in which the randomness enters linearly into the state equation as a force term, and a different optimization problem whose goal is to find a stochastic control u(ω) which has a similar sparsity pattern regardless of the realization ω. Notice further that the assumption ν > 0 does not eliminate the nonsmoothness of the objective functional, but it regularizes the optimal solution u, and is needed to use the fast optimation algorithm described in the following.
The well-posedness of (4.1) follows directly from standard variational arguments [39,17], being U ad a convex set, ϕ(u) := β u L 1 (D) a convex function and the objective functional coercive. In particular, the optimal solution u satisfies the variational inequality ([14, Proposition 2.2]) Through a pointwise discussion of the box constraints and an analysis of a Lagrange multiplier belonging to the subdifferential of ϕ in u, [37] showed that (4. where T : . Notice that F is nonsmooth due to the presence of the Lipschitz functions max(·) and min(·). Nevertheless, F can be shown to be semismooth [17], provided that T is continuously Fréchet differentiable, and further Lipschitz continuous interpreted as map from T : L 2 (D) → L r (D), with r > 2 [40,17]. These conditions are satisfied also in our settings since T is affine and further the adjoint variable p ω , solution of (2.5) with z = y d − S(u + f ), lies in , where r > 2 follows from the Sobolev embeddings.
Hence, to solve (4.3) we use the semismooth Newton method whose iteration reads for k = 1, 2, . . . until convergence, being the generalized derivative of F . Using the linearity of T and considering the supports of the weak derivatives of max(0, x) and min(0, x), we obtain that where χ is the charateristic function of the union of the disjoint sets It is possible to show that the generalized derivative G(u) is invertible with bounded inverse for all u, the proof being identical to the deterministic case treated in [36]. This further implies that the semismooth Newton method (4.4) converges locally superlinearly [40]. We briefly summarize these results in the following proposition.
Proposition 4.1. Let the inizialization u 0 be sufficiently close to the solution u of (4.1). Then the iterates u k generated by (4.4) converge superlinearly to u ∈ L 2 (D).
Introducing the supporting variables dy k ω and dp k w in L 2 (Ω; H 1 0 (D)), the semismooth Newton equation G(u k )du k = −F (u k ) may be rewritten as the equivalent saddle point system Further, if we set y 0 = S(f + u 0 ) and p 0 = S ⋆ (y d − y 0 ), due to the linearity of S and S ⋆ , it holds y k+1 = S(u k+1 ) = y k + dy k and similarly p k+1 = p k + dp k . Once fully discretized and using the notation E [p ω ] = N j=1 w j p j , the optimality condition (4.3) can be expressed through the nonlinear finite-dimensional map F : where the max(·) and min(·) functions act componentwise. Equation (4.6) leads to the saddle point system where H k ∈ R N h ×N h is a diagonal matrix representing the charateristic function To derive the expression of H, we assumed that a Lagrangian basis is used for the finite element space. Notice that (4.8) fits into the general form (1.2), and thus we Algorithm 4.1 Globalized semismooth Newton Algorithm to solve F(u) = 0 Require: u 0 , Tol ∈ R + , σ, ρ ∈ (0, 1).
use the collective multigrid algorithm to solve it. Further, with the notation of (1.2), it holds The collective multigrid iteration is then well-defined. The overall semismooth Newton Algorithm is summarized in Algorithm 4.1. At each iteration we solve (4.8) using the collective multigrid algorithm (line 3) and update the active sets given the new iteration (line 9). Notice that in order to globalize the convergence, we consider a line-search step (lines 5-7) performed on the merit function φ(u) = F(u) ⊤ M s F(u) [27].

Numerical experiments.
In this section we test the semismooth Newton algorithm for the solution of (4.3) and the collective multigrid algorithm to solve the related optimality system (4.8). We consider the random PDE-constraint (3.4) with the random diffusion coefficient (3.5). The semismooth iteration is stopped when φ(u k ) < 10 −9 . The inner linear solvers are stopped when the relative (unpreconditioned) residual is smaller than 10 −11 . Table 3 reports the number of semismooth Newton iterations and in brackets the averaged number of iterations of the V-cycle algorithm used as a solver (left) or as preconditioner for GMRES (right). Table 3 confirms the effectiveness of the multigrid algorithm, which requires essentially the same computational effort than in the linearquadratic case. Further, except for large values of σ 2 , we numerically observed that the line-search could be omitted without compromising the convergence of the outer iteration.
More challenging is the limit ν → 0 reported in Table 4. The performance of both the (globalized) semismooth Newton iteration and the inner multigrid solver Table 3 Number of semismooth Newton iterations, and average number of V-cycle iterations and preconditioned GMRES iterations (in brackets).   Table 4 Number of semismooth Newton iterations, and of V-cycle iterations and preconditioned GMRES iterations (in brackets). In the second row, the semismooth Newton method starts from a warm-up initial guess obtained through continuation. deteriorates. The convergence of the outer nonlinear algorithm can be improved by performing a continuation method, namely we consider a sequence of ν = 10 −j , j = 2, . . . , 8, and we start the j-th problem using as initial condition the optimal solution computed for ν = 10 −j+1 . Concerning the inner solver, the stand-alone multigrid algorithm struggles since for small values of ν the optimal control is of bang-bang type, that is satisfies u = a, u = b or u = 0 for almost every point of the mesh (for ν = 10 −8 only seven nodes are nonactive at the optimum). The matrices H k are then close to zero, and the multigrid hierarchy struggles to capture changes at such small scale. Nevertheless, the multigrid algorithm remains a very efficient preconditioner for GMRES even in this challenging limit. Fig. 1 shows a sequence of optimal controls for different values of β with and without box-constraints. The optimal control for β = 0 and without box-constraints corresponds to the minimizer of the linear-quadratic OCP (2.2). We observe that L 1 penalization indeed induces sparsity, since the optimal controls are more and more localized as β increases. Numerically we have verified that for sufficiently large β, the optimal control is identically equal to zero, a property shown in [37].

5.
A risk-adverse optimal control problem under uncertainty. In this section we consider an instance of risk-adverse OCPUU. This class of problems has recently drawn lot of attention since in engineering applications it is important to compute a control that minimizes the quantity of interest even in rare, but often troublesome, scenarios [20,21,1]. As a risk-measure [35], we use the Conditional Value-At-Risk (CVaR) of confidence level λ, λ ∈ (0, 1), that is, the expected value of a quantity of interest X given that the latter is greater than or equal to its λ-quantile, here denoted by VaR λ (X). Rockafellar and Uryasev [33] proved that CVaR λ (X) admits the equivalent formulation where (·) + := max(0, ·), if the distribution of X does not have an atom at VaR λ (X). In order to use tools from smooth optimization, we rely on a smoothing approach proposed in [20], which consists in replacing (·) + with a smooth function g ε , ε ∈ R + , such that g ε → (·) + in some functional norm as ε → 0. Specifically, we choose the C 2 -differentiable approximation Then, the smoothed risk-adverse OCPUU is where y ω ∈ V solves a ω (y ω , v) = (u + f, v) ∀v ∈ V, P-a.e. ω ∈ Ω.
where ν ∈ R + and λ ∈ [0, 1). The well-posedness of (5.3), the differentiability of its objective functional, as well as bounds for the error introduced by replacing (·) + with g ε (·), have been analyzed in [20]. Further, defining Q ω = 1 2 y ω − y d 2 L 2 (D) − t, the optimality conditions form the nonlinear system, Approximating V and E with V h and E, and letting x = (y, u, p, t), the finitedimensional discretization of (5.4) correponds to the nonlinear system for j = 1, . . . , N .
A possible approach to solve (5.5) is to use a Newton method, which given (5.9) for i = 1, . . . , N . Unfortunately, J k can be singular away from the optimum, in particular whenever E g ′′ ε (Q k ω ) = 0 which implies . . , N, (5.10) which is not unlikely for small ε since supp(g ′′ ε ) = (− ε 2 , ε 2 ). Splitting strategies have been proposed (e.g. [26] in a reduced approach), in which whenever (5.10) is satisfied, an intermediate value of t is computed by solving F 4 (t; y, u, p) = 0 so to violate (5.10). In the next section, we discuss a similar splitting approach. To speed up the convergence of the outer nonlinear algorithm, we use a preconditioned Newton method based on nonlinear elimination [22]. At each iteration we will need to invert saddle-point matrices like (1.2), possibly several times. To do so, we rely on the collective multigrid algorithm.

Nonlinear preconditioned Newton method.
Nonlinear elimination is a nonlinear preconditioning technique based on the identification of variables and equations of F (e.g. strong nonlinearities) that slow down the convergence of Newton method. These components are then eliminated through the solution of a local nonlinear problem at every step of an outer Newton. This elimination step provides a better initial guess for the outer iteration, so that a faster convergence is achieved [22,41].
In light of the possible singularity of J, we split the discretized variables x into x = (x, t), and we aim to eliminate the variables x to obtain a scalar nonlinear equation only for t. To do so, we partition (5.4) as )N h , and J 2,2 ∈ R. Notice that J 1,1 is always nonsingular, while J 2,1 , J 1,2 and J 2,2 are identically zero if (5.10) is verified. Thus F 1 allows us to define an implicit map h : R → R (2N +1)N h , such that F 1 (h(t), t) = 0, so that the first set of nonlinear equations in (5.11) are satisfied. We are then left to solve the nonlinear scalar equation To do so using the Newton method, we need the derivative of F (t) evaluated at t = t k which, using implicit differentiation, can be computed as The nonlinear preconditioned Newton method is described in Alg. 5.1, and consists in solving (5.12) with Newton method. However, to overcome the possible singularity of J k 2,2 , J k 1,2 and J k 2,1 , we check at each iteration k if (5.10) is satisfied, and in the affermative case we update x k by solving F 1 (x k+1 , t k ) = 0 using Newton method, and update t k by solving F 2 (x k , t k+1 ) = 0. Notice further, that each iteration of the backtracking line-search requires to solve F 1 (h(t), t) = 0 using Newton method, thus additional linear systems with matrix J 1,1 must be solved.
We report that we also tried to eliminate t by computing the map l such that F 2 (x, l(x)) = 0, while iterating on the variable x. This has the advantage that l can be evaluted very cheaply, being a scalar equation. However, we needed many more Algorithm 5.1 Nonlinear preconditioned Newton method to solve F( x) = 0.
iterations both of the outer Newton method, and consequently of the inner linear solver. Thus, according to our experience, this second approach resulted less efficient and appealing.

Numerical experiments.
In this section we report numerical tests to asses the performance of the preconditioned Newton algorithm to solve (5.5), and of the collective multigrid algorithm to invert the matrix J 1,1 . We consider the random PDEconstraint (3.4) with the random diffusion coefficient (3.5). Table 5 reports the number of outer and inner Newton iterations, and the average number of V-cycle iterations and of preconditioned GMRES iterations to solve the linear systems at each (inner/outer) Newton iterations. The outer Newton iteration is stopped when |F (t k )| ≤ 10 −6 , the inner Newton method to compute h(·) is stopped when F 1,1 (x; t k ) 2 ≤ 10 −9 , and the linear solvers are stopped when the (unpreconditioned) residual is smaller than 10 −10 .
In Table 5, the number of outer Newton iterations is quite stable, while the number of inner Newton iterations varies between five and twenty iterations per outer iteration, essentially due to how difficult is to compute the nonlinear map h(t) by solving F 1 (x; t k ) = 0 in line (5), (8) and (11) of Alg. 5.1. The average number of inner linear solvers is quite stable as well. We emphasize that the top left blocks of J 1,1 involve the matrices C i (y k i , t k ) (see (5.8)) which contain a dense low-rank term if g ′′ ε (Q k ωi ) = 0. Further, as ε → 0, g ′′ ε (·) tends to a Dirac delta, so the dense term become dominant. Multigrid methods based on a pointwise relaxations are expected to be not very efficient for these matrices which may not be diagonally dominant. The standard V-cycle algorithm indeed struggles, however the Krylov acceleration performs much better as it handles these low-rank perturbation with smaller effort. We sometimes noticed that the GMRES residual stagnates after 20/30 iterations around 10 −8 /10 −9 , due to a loss of orthogonality in the Krylov subspace. We allowed a maximum number of 80 GMRES iterations per linear system. Figure 2 compares the two optimal controls obtained minimizing either E [Q(y ω )] or CVaR λ (Q(y ω )), and the quantity of interests Q(y ωj ) for every sample discretization Table 5 Number of preconditioned Newton iterations, and of V-cycle iterations and preconditioned GM-RES iterations (in brackets). at optimum. We remark that the risk-adverse control indeed minimizes the risk of having large values of Q(y ω ), at the price of leading to higher values of Q(y ω ) for most realizations. Notice that we used Gauss-Legendre nodes, so that not all realizations have the same weights. The more likely events are around the median, while the first 17 realizations count for around 1% of the discretized probability mass.
6. Conclusion. We have presented a multigrid method to solve the large saddle point linear systems that typically arise in full-space approaches to solve OCPUU.
We tested the algorithm as a iterative solver and as a preconditioner on three test cases: a linear-quadratic OCPUU, a nonsmooth OCPUU, and a risk-adverse nonlinear OCPUU. Overall, the multigrid method shows very good performances and robustness with respect to the several parameters of the problems considered.