Optimal design with uncertainties: a risk-averse approach

Optimal design with uncertain ties: a risk-a v erse approac h Amal Alphonse ∗ P etar Kun ˇ stek † Mark o V rdoljak ‡ F ebruary 24, 2026 Abstract W e study a class of sto c hastic optimal design problems for elliptic partial differen tial equations in div ergence form, where the coefficients represen t mixtures of tw o conducting materials. The ob jectiv e is to minimize a generalized risk measure of the system resp onse, incor- p orating uncertaint y in the loading through probabilit y distributions. W e establish existence of relaxed optimal designs via homogenization theory and deriv e first-order stationarity conditions satisfied by the optima. Based on these necessary conditions, we develop an opti- malit y criteria algorithm for n umerical computations. The sto chastic comp onen t is treated using a truncated Karh unen–Lo` ev e expansion, allo wing ev aluation of the v alue-at-risk (V aR) and conditional v alue- at-risk (CV aR) con tributions arising from the sensitivit y analysis and featured in the algorithm. The metho d is illustrated for an example in volving CV aR-based compliance minimization. 1 In tro duction Optimal design of structures and materials gov erned b y partial differen tial equations has b een extensively studied in the deterministic setting, where uncertain ties in loading, b oundary conditions, and material prop erties are neglected [ 22 , 6 , 1 , 5 , 25 ]. In many applications, ho wev er, such uncertainties ∗ W eierstrass Institute, Berlin, Germany , alphonse@wias-berlin.de † Departmen t of Mathematics, F aculty of Science, Univ ersity of Zagreb, Croatia, petar@math.hr ‡ Departmen t of Mathematics, F aculty of Science, Univ ersity of Zagreb, Croatia, marko@math.hr 1 ha ve a significan t influence on the system response, and deterministic opti- mization ma y lead to designs that p erform p o orly or ev en fail under random p erturbations. This has partly motiv ated the developmen t of sto chastic and risk-a verse PDE-constrained optimization, in which random inputs are ex- plicitly mo deled and optimization criteria are form ulated to control statistical prop erties of the resulting random state [ 10 , 2 , 11 , 21 ]. A foundational framework for risk-av erse PDE-constrained optimization w as established b y Kouri and Suro wiec in [ 16 , 18 ] (see also [ 17 ]). They studied optimization problems constrained b y elliptic PDEs with random forcin g and incorp orated coheren t risk measures, most notably the Conditional V alue-at- Risk (CV aR) to quantify the risk of high-impact, lo w-probability outcomes. Their analysis included the existence of optimal con trols and the deriv ation of first-order optimalit y conditions. The presen t w ork is motiv ated by extending these ideas to optimal de- sign problems arising in comp osite materials, where the con trol v ariable is the spatial distribution of t w o phases with differen t conductivity prop erties. Incorp orating uncertain ty into such design problems presen ts new mathe- matical c hallenges. The resp onse of the system to random loads b ecomes a random field, and the ob jectiv e is t ypically a functional of the en tire dis- tribution, not merely its exp ectation. As highlighted in e.g., [ 24 , 14 ] and references therein, coherent risk measures (suc h as CV aR) pro vide a robust alternativ e to exp ected-v alue ob jectives b y p enalizing rare but severe realiza- tions. How ev er, the presence of risk measures induces nonsmo othness in the ob jective functional, which in turn requires careful mathematical analysis. Let us presen t the problem of interest in this w ork. The aim is to find the b est arrangemen t of giv en materials in an op en and bounded domain Ω ⊆ R d so that the resulting b o dy has optimal prop erties regarding m different regimes. W e consider the simplest case of tw o isotropic materials p ossessing conductivities α and β with 0 < α < β . If we let Ω α ⊆ Ω be the measurable set occupied by the first material, and χ ≡ χ Ω α its c haracteristic function, then the conductivit y can b e written as A = χα I + (1 − χ ) β I (1) where I is the identit y . Let ( S, F , P ) b e a probabilit y space, and supp ose that f i ∈ L p ( S ; H − 1 (Ω)) are giv en for some p ∈ (1 , ∞ ), and introduce the b oundary v alue problems (on a p oin twise a.s. lev el)  − div( A ∇ u i ) = f i ( s, · ) u i ∈ H 1 0 (Ω) i = 1 , . . . m . (2) F or example, u i could mo del the electrical p oten tial or temperature for the i th state. In general, these equations allow us to take into account m different 2 states and reflect the fact that a device is often required to b e robust in mul- tiple scenarios. Define a quantit y of interest b J : L ∞ (Ω; { 0 , 1 } ) → L 0 ( S ; R ) (where b y L 0 w e mean the set of measurable functions) that measures the efficacy of the mixture by b J ( χ )( s ) = Z Ω  χ ( x ) g α ( s, x , u ( s, x )) + (1 − χ ( x )) g β ( s, x , u ( s, x ))  d x where u = ( u 1 , . . . , u m ) with the u i determined by ( 2 ), A is giv en b y ( 1 ) and g α and g β are sufficiently smo oth functions satisfying certain gro wth conditions, see Theorem 1.1 for some examples. Although g α and g β could dep end on u i in a more complicated fashion (such as dep ending on ∇ u i ) and suc h cases are also interesting in applications, in this work, w e consider only dep endence on the p oin t wise state v alues. The functional is considered in this form in [ 22 ] for the single-state prob- lem and in [ 1 ] for the m ulti-state case, as it encompasses the most general class under consideration, while also allo wing for a p ossible dep endence of the functional on the random parameter s . In order to treat the fact that the ob jective function is probabilistic, w e can consider a risk-a verse optimal design problem. F or a giv en risk mea- sure R : L 1 ( S ) → R , the aim is to find a characteristic function χ on Ω that minimizes b J in a risk-av erse sense: min χ ∈ L ∞ (Ω; { 0 , 1 } ) R [ b J ( χ )], p ossibly supplemen ted by some regularization term. As is standard in deterministic optimal design, in order to a void ill p osed- ness of the problem (see e.g. [ 1 , § 3.1.5]), a prop er relaxation of the problem is required, leading to the use of comp osite materials formed as fine mix- tures of the original phases. Accordingly , instead of characteristic functions for the first material, its lo cal v olume fraction θ , taking v alues in [0 , 1], is in tro duced. Ho wev er, the lo cal volume fraction θ alone does not fully define the comp osite; rather a set of admissible comp osites whic h we denote by K ( θ ) (this contains all the effective materials with prop ortions θ ) do es, see ( 8 ). Therefore, the con trol v ariable becomes ( θ , A ) with A ∈ K ( θ ) almost ev erywhere in Ω. Let us denote A := n ( θ , A ) ∈ L ∞ (Ω; [0 , 1]) × L ∞ (Ω; Sym d ) : A ( x ) ∈ K ( θ ( x )) for a.e. x ∈ Ω o (3) where the notation Sym d means the set of real symmetric d × d matrices. Finally , the problem w e study in this pap er can b e form ulated as min ( θ, A ) ∈A R [ J ( θ , A )] +  ( θ ) (P) 3 where J ( θ , A )( s ) = Z Ω  θ ( x ) g α ( s, x , u ( s, x ))+(1 − θ ( x )) g β ( s, x , u ( s, x ))  d x a.s. in s, (4) u is uniquely determined b y A through ( 2 ), R is a risk measure and  is a regularizer. W e will form ulate our problem in full rigor in Section 2 . Let us giv e some examples of the ob jectiv e functional, risk measure and regularizer. Example 1.1. A n example of the obje ctive function for m = 1 is obtaine d with the choic e g α ( s, x , u ) = g β ( s, x , u ) = f ( s, x ) u , whenc e J ( θ , A )( s ) = Z Ω f ( s, x ) u ( s, x ) d x , which is the ener gy dissip ation functional. Another example r elates to tr acking- typ e obje ctives typic al in optimal c ontr ol wher e g α ( s, x , u ) = | u ( s, x ) − u α ( x ) | 2 and g β ( s, x , u ) = | u ( s, x ) − u β ( x ) | 2 wher e u α , u β ar e given (deterministic) de- sir e d states; this choic e le ads to J ( θ , A )( s ) = Z Ω ( θ ( x ) | u ( s, x ) − u α ( x ) | 2 + (1 − θ ( x )) | u ( s, x ) − u β ( x ) | 2 d x . F or the multiple state pr oblem with m > 1 , we c ould c onsider simply the sum of individual functions asso ciate d to e ach state: J ( θ , A )( s ) = m X i =1 Z Ω  θ ( x ) g i α ( s, x , u i ( s, x )) + (1 − θ ( x )) g i β ( s, x , u i ( s, x ))  d x . These examples wer e p ointe d out in [ 1 , § 3.1]. R e gar ding the risk me asur e, we c an take R to b e CV aR as mentione d ab ove (se e ( 23 ) for a definition) or we c ould work in the risk-neutr al setting wher e R = E is the exp e cte d value, for example. F or the r e gularizer  , one c ould c onsider  ( θ ) = 1 2 ∥ θ ∥ 2 L 2 (Ω) or the typic al choic e  ( θ ) = λ R Ω θ d x . Remark 1.2. A differ ent appr o ach is to p ose the pr ob ability maximization pr oblem max ( θ, A ) ∈A P ( J ( θ , A )( s ) ≤ c ) (5) for a given ac c eptable thr eshold level c > 0 . This mo dels the situation wher e we se ek to find the design that maximizes the pr ob ability that our obje ctive is satisfie d up to a smal l err or. This is of c ourse a differ ent mo del, but it is not to o difficult to se e that it c an b e written in the form ( P ) with the risk 4 me asur e R chosen to b e V aR γ . However, our the ory do es not apply her e b e c ause V aR γ is not c onvex and c onvexity is enfor c e d on R in The or em 2.4 . Nevertheless, pr ob ability maximization pr oblems such as ( 5 ) c an b e tackle d dir e ctly and stationarity c onditions c an b e derive d, at le ast in the setting wher e the optimization is unc onstr aine d and some structur al assumptions on the obje ctive in terms of the unc ertain variable ar e pr esent (such as c onvexity or quadr aticity, the latter of which is satisfie d by our example in Se ction 4 ). Differ entiability of the pr ob ability function has b e en the obje ct of study in e.g. [ 27 , 29 ] using the spheric al r adial de c omp osition, mostly in the finite- dimensional setting and to a lesser de gr e e the r eflexive, sep ar able Banach sp ac e setting [ 13 ], which do es not fit our non-sep ar able, non-r eflexive c ontr ol sp ac e A . The work [ 28 ] do es c onsider a gener al Banach sp ac e setting, but utilizes the ab ove-mentione d c onvexity. We b elieve it is p ossible to extend these r esults to our c onstr aine d optimization pr oblem in a non-r eflexive, non- sep ar able sp ac e with a nonc onvex dep endenc e on the unc ertainty (sinc e the pr o ofs in the ab ove-cite d p ap ers either c an b e adapte d or they simply hold verb atim when the obje ctive function is C 1 , as it is for us), yet this is non- trivial and outside the sc op e of this p ap er. Thus we le ave this, and other inter esting asp e cts such as p ointwise chanc e c onstr aints, for investigation in futur e work. W e conclude this section with a brief ov erview of the pap er. In Sec- tion 2 , we precisely introduce the set of admissible controls and state the assumptions that ensure the existence of a solution. W e also outline the k ey elemen ts of the homogenization theory employ ed. Section 3 is devoted to the deriv ation of first-order optimality conditions. Finally , Section 4 presents a v arian t of the optimality criteria algorithm tailored to the n umerical solution of the problem. By emplo ying the Karhunen–Lo ` ev e expansion of the random righ t-hand sides, w e obtain an explicit finite-dimensional representation of the sto c hastic component, whic h enables an efficien t implementation of the algorithm and is illustrated by a n umerical example. W e conclude the pap er in Section 5 with some remarks. 2 H -con v ergence and existence of solutions Here and b elo w, w e embed the Sob olev space H 1 0 (Ω) with the norm ∥ u ∥ H 1 0 (Ω) := ∥∇ u ∥ L 2 (Ω) . In connection with the problem inv olving b J (wherein w e note that χ influences A from ( 1 ), which in turn influences the b oundary v alue problems ( 2 )), it is well known that no reasonable pair of topologies exists for which 5 the mapping χ 7→ u is con tinuous from the set of measurable characteristic functions on Ω to H 1 (Ω), where u denotes the unique solution of  − div( A ∇ u ) = f u ∈ H 1 0 (Ω) , (6) with A given b y ( 1 ). On the other hand, homogenization theory [ 23 , 22 , 1 ] pro vides a natural topology on a class of co efficients A for whic h the mapping A 7→ u is con tinuous for every f ∈ H − 1 (Ω) with u equipp ed with the w eak H 1 top ology . T o b e more precise, with M d ( R ) denoting the space of real d × d square matrices, this topology is introduced on the set M ( α, β ; Ω), defined [ 23 , 22 ] M ( α, β ; Ω) := { A ∈ L ∞ (Ω; M d ( R )) : A ( x ) ξ · ξ ≥ α | ξ | 2 and A ( x ) ξ · ξ ≥ 1 β | A ( x ) ξ | 2 for a.e. x ∈ Ω and for all ξ ∈ R d } , where | · | denotes the Euclidean norm on R d . The notion of H -con v ergence is imp ortan t, whic h w e recall now. Definition 2.1 ( H -con vergence) . We say that the matrix se quenc e { A n } ⊆ M ( α, β ; Ω) H -con verges to A ⊆ M ( α , β ; Ω) with r esp e ct to the H -top olo gy if for arbitr ary f ∈ H − 1 (Ω) the solutions u n ∈ H 1 0 (Ω) of − div( A n ∇ u n ) = f satisfy u n  u in H 1 0 (Ω) and A n ∇ u n  A ∇ u in L 2 (Ω; R d ) wher e u ∈ H 1 0 (Ω) solves − div( A ∇ u ) = f . We write this as A n H  A . Crucially , this con vergence in tro duces a top ology that is kno wn to b e metrizable and compact on M ( α, β ; Ω) [ 26 , Theorem 6]. Note also that the set of conductivities of the form ( 1 ) where χ is an arbitrary c haracteristic function b elongs to M ( α, β ; Ω). If a sequence of c haracteristic functions χ n on Ω w eakly- ∗ conv erges to θ ∈ L ∞ (Ω; [0 , 1]) and A n = χ n α I + (1 − χ n ) β I H -con verges to A ∈ M ( α, β ; Ω), w e sa y that A is a homo genize d c onductivity of a t wo-phase comp osite material obtained by mixing α I and β I in prop ortions θ and 1 − θ , see [ 1 , § 2.1.2]. The set of all suc h homogenized conductivities A (i.e., H -con vergen t limits A of the pro cess describ ed ab o ve) is denoted by K θ : K θ = { A ∈ M ( α, β ; Ω) : A is the H -limit of A n as describ ed ab ov e . } (7) 6 This set can b e understo o d as the set of all comp osite materials associated to θ . By the lo cal character of H -conv ergence it follo ws that the set of homog- enized conductivities K θ can be c haracterized as L ∞ matrix functions such that for almost ev ery x ∈ Ω w e ha ve A ( x ) ∈ K ( θ ( x )). Here, for an y ϑ ∈ [0 , 1] the set K ( ϑ ) is some closed set of matrices. In fact, in our situation of mix- ing t wo isotropic phases, K ( ϑ ) is a set of symmetric matrices which can b e explicitly characterized in terms of their eigenv alues (see [ 22 , Prop osition 10] for precise details, see also [ 20 ], [ 1 , Theorem 2.2.13]): K ( ϑ ) = ( M ∈ Sym d : λ − ϑ ≤ λ i ≤ λ + ϑ for ev ery 1 ≤ i ≤ d, d X i =1 1 λ i − α ≤ 1 λ − ϑ − α + d − 1 λ + ϑ − α , d X i =1 1 β − λ i ≤ 1 β − λ − ϑ + d − 1 β − λ + ϑ ) (8) where w e used the notation λ i to mean the eigenv alues of the matrix M and λ − ϑ and λ + ϑ are the harmonic and arithmetic means of α and β , i.e., λ − ϑ = ( ϑ/ α + (1 − ϑ ) /β ) − 1 and λ + ϑ = ϑα + (1 − ϑ ) β . No w it is easy to see that the set A 1 defined in ( 3 ) A := n ( θ , A ) ∈ L ∞ (Ω; [0 , 1]) × L ∞ (Ω; Sym d ) : A ( x ) ∈ K ( θ ( x )) for a.e. x ∈ Ω o is compact with resp ect to the weak- ∗ top ology for θ and the H -top ology for A . F or fixed ϑ, the set K ( ϑ ) is con v ex [ 1 , Theorem 2.2.13], but the set A is not [ 32 , Example 2.1]. W e return to the problem ( P ). T o b egin with, let us consider the well p osedness of the constituent state equations ( 2 ). Lemma 2.2. F or A ∈ M ( α, β ; Ω) , the pr oblem ( 2 ) is uniquely solvable and u i ∈ L p ( S ; H 1 0 (Ω)) for i = 1 , . . . , m . F urthermor e, we have the fol lowing a.s. estimate: ∥ u i ( s, · ) ∥ H 1 0 (Ω) ≤ α − 1 ∥ f i ( s, · ) ∥ H − 1 (Ω) , i = 1 , . . . , m . (9) Pr o of. W ell p osedness almost surely in s ∈ S , together with the a priori estimate ( 9 ), follows directly from the Lax–Milgram lemma. Defining the linear op erator T : H 1 0 (Ω) → H − 1 (Ω) by T u := − div ( A ∇ u ), w e can write u i ( s ) = T − 1 f i ( s ). By [ 34 , Section V.5, Corollary 2] u i is strongly measurable and, b y the estimate ( 9 ) ab o ve, b elongs to the sp ecified Bo chner space. 1 Note that A corresp onds to the set C D defined in [ 1 , Equation (3.39)]. 7 Since the righ t-hand sides f 1 , . . . , f m ∈ L p ( S ; H − 1 (Ω)) from ( 2 ) are con- sidered fixed, to simplify notation we introduce the mapping G : M ( α , β ; Ω) →  L p ( S ; H 1 0 (Ω))  m , G ( A ) = u determined b y ( 2 ). Recall that the (relaxed) function J : A → L 1 ( S ) from ( 4 ) is defined via J ( θ , A )( s ) = Z Ω  θ ( x ) g α ( s, x , u ( s, x )) + (1 − θ ( x )) g β ( s, x , u ( s, x ))  d x suc h that u = G ( A ), i.e. u and A are related via ( 2 ). Regarding the functions g α , g β , w e assume the follo wing. Assumption 2.3. The functions g α and g β ar e Car ath ´ eo dory functions (me a- sur able in ( s, x ) on the pr o duct sp ac e S × Ω and c ontinuous in u ) satisfying the gr owth c onditions: | g γ ( s, x , u ) | ≤ ϕ γ ( s, x ) + ψ γ ( s, x ) | u | q 1 for γ = α, β , (10) with 1 ≤ q 1 ≤ min { p, q ∗ } wher e q ∗ =  + ∞ : d ≤ 2 , 2 d d − 2 : d > 2 , and ϕ γ ∈ L 1 ( S × Ω) , ψ γ ∈ L q 2 ( S ; L q 3 (Ω)) with q 2 = p p − q 1 , and q 3 > q ∗ q ∗ − q 1 . In addition, we make the following assumption 2 on the risk measure R and the regularization term  . Assumption 2.4. L et R : L 1 ( S ) → R b e lower semic ontinuous and c onvex and let  : L ∞ (Ω; [0 , 1]) → R b e we ak- ∗ lower semic ontinuous. Under these assumptions, the relaxed problem ( P ), which w e recall again min ( θ, A ) ∈A R [ J ( θ , A )] +  ( θ ) , ( P ) is a w ell defined problem. Let us now show that it p ossesses a solution. Theorem 2.5. Supp ose that The or ems 2.3 and 2.4 ar e satisfie d. Then, the r elaxe d pr oblem ( P ) admits a solution ( θ ∗ , A ∗ ) ∈ A . 2 Actually , the conv exity assumption on R is only used from Theorem 3.5 onw ards. 8 Pr o of. Let us first sho w that R ◦ J : A → R is sequentially low er semicon- tin uous. Consider a sequence ( θ n , A n ) ∈ A con verging to ( θ , A ) with resp ect to the appropriate top ologies, so that θ n ∗  θ in L ∞ (Ω) and A n H  A . Let us set u n = G ( A n ) and u = G ( A ). By H -con vergence we hav e the weak con vergence u n i ( s, · )  u i ( s, · ) in H 1 0 (Ω) , a . s . (11) Moreo ver, b y the Sob olev imbedding theorem we hav e H 1 0 (Ω) ⊂ L q (Ω) com- pactly for any q ∈ [1 , q ∗ ) with q ∗ as defined ab ov e. Hence the conv ergence is actually strong in L q (Ω). F urthermore, a.s. in s ∈ S , applying the a priori estimate ( 9 ) for ( 2 ), we ha ve ∥ u n i ( s, · ) ∥ H 1 0 (Ω) ≤ α − 1 ∥ f i ( s, · ) ∥ H − 1 (Ω) , i = 1 , . . . , m . (12) Using the Sobolev embedding mentioned ab o v e, we further obtain a.s. ∥ u n i ( s, · ) ∥ L q (Ω) ≤ C α − 1 ∥ f i ( s, · ) ∥ H − 1 (Ω) , i = 1 , . . . , m for any q ∈ [1 , q ∗ ). The right-hand side of this inequalit y (for an y i ) b e- longs to L p ( S ). Combined with the strong L q (Ω) conv ergence of u n i ( s, · ) (for q ∈ [1 , q ∗ )) and Leb esgue’s dominated con vergence theorem [ 15 , Prop osition 1.2.5], this implies the strong con vergence of u n to u in L p ( S ; L q (Ω; R m )). Consequen tly , there exists a function U ∈ L p ( S ; L q (Ω)) suc h that, for a sub- sequence that w e shall relab el, | u n | ≤ U almost surely-almost ev erywhere on S × Ω. It follo ws that u n → u in L min( p,q ) ( S × Ω; R m ). Hence, up to a further subsequence, u n ( s, x ) → u ( s, x ) for a.s.-a.e. ( s, x ) ∈ S × Ω, and therefore g γ ( s, x , u n ( s, x )) → g γ ( s, x , u ( s, x )) , γ ∈ { α, β } , a.s.-a.e. ( s, x ) ∈ S × Ω , (13) since g α and g β are Carath ´ eo dory functions. By the gro wth conditions ( 10 ), w e hav e | g γ ( s, x , u n ( s, x )) | ≤ ϕ γ ( s, x ) + ψ γ ( s, x ) U ( s, x ) q 1 . (14) Since U q 1 ∈ L p q 1 ( S ; L q q 1 (Ω)), an application of H¨ older’s inequality to the second term on the right-hand side (see [ 12 , Lemma IV.1.4], or the analogus result for the mixed-norm Leb esgue spaces [ 7 ]) yields ψ γ U q 1 ∈ L 1 ( S × Ω) if 1 q 2 + q 1 p = 1 , 1 q 3 + q 1 q = 1 whic h hold by assumption. Since the right-hand side of ( 14 ) b elongs to L 1 ( S × Ω), the p oin twise con- v ergence in ( 13 ), together with the dominated con v ergence theorem, implies that g γ ( · , · , u n ) con verges to g γ ( · , · , u ) strongly in L 1 ( S × Ω), for γ ∈ { α , β } . 9 Finally , since θ n ∗  θ in L ∞ (Ω) we conclude J ( θ n , A n ) → J ( θ , A ) in L 1 ( S ). Indeed, we observe that Z S |J ( θ n , A n ) − J ( θ, A ) | d P ( s ) = Z S     Z Ω  θ n g α ( u n ) + (1 − θ n ) g β ( u n )  −  θ g α ( u ) + (1 − θ ) g β ( u )  d x     d P ( s ) = Z S    Z Ω  θ n ( g α ( u n ) − g α ( u )) + (1 − θ n ) ( g β ( u n ) − g β ( u ))  d x + Z Ω  ( θ n − θ ) g α ( u ) + ( θ − θ n ) g β ( u )  d x    d P ( s ) ≤ Z S Z Ω | θ n ( g α ( u n ) − g α ( u )) | d x d P ( s ) + Z S Z Ω | ( θ n − θ ) g α ( u ) | d x d P ( s ) + Z S Z Ω | (1 − θ n ) ( g β ( u n ) − g β ( u )) | d x d P ( s ) + Z S Z Ω | ( θ − θ n ) g β ( u ) | d x d P ( s ) = I + I I + I I I + IV . No w, we hav e I = Z S Z Ω | θ n ( g α ( u n ) − g α ( u )) | d x d P ( s ) ≤ ∥ θ n ∥ L ∞ (Ω) ∥ g α ( u n ) − g α ( u ) ∥ L 1 ( S × Ω) → 0 since ( θ n ) is uniformly b ounded in L ∞ due to w eak- ∗ conv ergence, and the second factor tends to zero by the strong L 1 ( S × Ω) con vergence established ab o ve. Next, consider the term I I ab o ve. Using the weak- ∗ conv ergence, we find h n ( s ) := Z Ω | ( θ n − θ ) g α ( u )( s ) | d x → 0 p oin twise a.s. in s . Moreo ver, we hav e the b ound h n ( s ) ≤ C Z Ω | g α ( u )( s ) | d x ≤ C Z Ω ϕ γ ( s ) + ψ γ ( s ) | u ( s ) | q 1 d x b y the growth condition ( 10 ). The right-hand side belongs to L 1 ( S ) by assumption, as explained previously . Therefore, the dominated conv ergence theorem implies that h n → 0 in L 1 ( S ), i.e., R S h n ( s ) → 0, which sho ws that I I conv erges to 0. The terms I I I and IV can b e treated similarly . This establishes that J ( θ n , A n ) → J ( θ , A ) in L 1 ( S ). 10 As men tioned b efore, the top ology introduced on A is metrizable, whic h enables us to conclude that J : A → L 1 ( S ) is con tin uous. It follo ws that the sum F := R ◦ J +  is also lo wer semicontin uous with respect to w eak −∗ top ology for θ and H -top ology for A . No w the result follows b y the W eier- strass theorem, see [ 19 , § 7, Theorem 7.3.1]. 3 Necessary conditions of optimalit y In this section w e will derive first-order stationarity conditions for ( P ). W e will denote by ( θ ∗ , A ∗ ) an optimal relaxed solution of ( P ), the existence of whic h is assured by Theorem 2.5 . W e need additional regularity on the functions g α and g β in the form of the next assumption. Assumption 3.1. The functions ∂ g α ∂ u i and ∂ g β ∂ u i ar e supp ose d to b e Car ath´ eo dory functions (me asur able in ( s, x ) and c ontinuous in u ) satisfying the gr owth c onditions     ∂ g γ ∂ u i ( s, x , u )     ≤ e ϕ γ ( s, x ) + e ψ γ ( s, x ) | u | q 1 − 1 for γ = α, β ; i = 1 , . . . , m, (15) wher e e ϕ γ ∈ L p ′ ( S ; L q ′ ∗ (Ω)) , e ψ γ ∈ L q 4 ( S ; L q 5 (Ω)) with 1 ≤ q 1 ≤ min { p, q ∗ } , q 4 = p p − q 1 , and q 5 = q ∗ q ∗ − q 1 . In [ 22 ] (and similarly [ 1 , Theorem 3.2.13] for multiple state problems), the Gˆ ateaux differen tial of the functional J ( · , · )( s ) : A → R is computed for a.s. s ∈ S . Using this and arguing lik e in the pro of of Theorem 2.5 , w e obtain the Gˆ ateaux differen tiability of J : A → L 1 ( S ) (see also [ 33 ]): J ′ ( θ ∗ , A ∗ )( δ θ , δ A )( s ) = Z Ω δ θ  g α ( s, x , u ∗ ( s, x )) − g β ( s, x , u ∗ ( s, x ))  d x − Z Ω m X i =1 δ A ∇ u ∗ i ( s, x ) · ∇ p ∗ i ( s, x ) d x , (16) where p ∗ = ( p ∗ 1 , . . . , p ∗ m ) denotes the adjoin t state, defined as the unique solution 3 of ( − div( A ∗ ∇ p i ) = θ ∗ ∂ g α ∂ u i ( · , · , u ∗ ) + (1 − θ ∗ ) ∂ g β ∂ u i ( · , · , u ∗ ) p i ∈ L p ′ ( S ; H 1 0 (Ω)) i = 1 , . . . , m. (17) 3 Existence and uniqueness for ( 17 ) can b e prov ed in a similar fashion as Theorem 2.5 ; w e omit the details. 11 In fact, w e will prov e that J is differentiable in the F r´ echet sense (see also [ 22 , Prop osition 7]). This will b e used later when we come to apply the c hain rule. F or this purp ose, we need the following technical lemma. Lemma 3.2. If A n → A in L ∞ (Ω; M α,β ) and f n → f ∈ H − 1 (Ω) , the se quenc e ( u n ) of solutions to  − div( A n ∇ u n ) = f n u n ∈ H 1 0 (Ω) c onver ges str ongly to u in H 1 0 (Ω) , wher e u is the solution of ( 6 ) . Pr o of. The b oundedness of Ω implies that A n → A in L 1 (Ω; M d ( R )), and hence, by [ 1 , Lemma 1.2.22], A n → A in the sense of H -conv ergence. Hence u n  u in H 1 0 (Ω). T o establish strong conv ergence, first write − div( A ∇ u n ) = f n − div (( A − A n ) ∇ u n ) =: g n . Since ∇ u n  ∇ u in L 2 (Ω), the sequence ( ∇ u n ) is b ounded in L 2 (Ω), and w e deduce that ( A − A n ) ∇ u n → 0 in L 2 (Ω; R d ). This implies ∇ · ( A − A n ) ∇ u n → 0 in H − 1 (Ω), i.e., g n → f in H − 1 (Ω). Finally , since the mapping u 7→ − div( A ∇ u ) is an isomorphism betw een H 1 0 (Ω) and H − 1 (Ω), w e conclude that u n → u strongly in H 1 0 (Ω). Lemma 3.3. L et The or ems 2.3 , 2.4 and 3.1 hold. Then, the mapping J : A → L 1 ( S ) is c ontinuously differ entiable. Pr o of. Let A n → A strongly in L ∞ (Ω; M α,β ) and θ n → θ strongly in L ∞ (Ω), and denote the corresp onding states u n = G ( A n ) and u = G ( A ). Our goal is to sho w that J ′ ( θ n , A n ) con v erges to J ′ ( θ , A ) strongly in the operator norm. 12 T o this end, observ e that ∥ ( J ′ ( θ n , A n ) − J ′ ( θ ∗ , A ∗ ))( δ θ , δ A ) ∥ L 1 ( S ) = Z S    Z Ω δ θ  g α ( u n ) − g β ( u n ) − ( g α ( u ∗ ) − g β ( u ∗ ))  d x − Z Ω m X i =1 δ A ∇ u n i · ∇ p n i − m X i =1 δ A ∇ u ∗ i · ∇ p ∗ i d x    d P ( s ) ≤ Z S Z Ω | δ θ  g α ( u n ) − g β ( u n ) − ( g α ( u ∗ ) − g β ( u ∗ ))  | d x d P ( s ) + Z S Z Ω m X i =1 | δ A ∇ u n i · ∇ p n i − δ A ∇ u ∗ i · ∇ p ∗ i | d x d P ( s ) ≤ ∥ δ θ ∥ L ∞ (Ω) Z S Z Ω | g α ( u n ) − g β ( u n ) − ( g α ( u ∗ ) − g β ( u ∗ )) | d x d P ( s ) + ∥ δ A ∥ L ∞ (Ω;Sym d ) Z S Z Ω m X i =1 |∇ u n i · ∇ p n i − ∇ u ∗ i · ∇ p ∗ i | d x d P ( s ) . The first integral on the righ t-hand side can b e shown to v anish by following the argumen t in the pro of of Theorem 2.5 , whic h yields the strong con ver- gence g α ( · , · , u n ) − g β ( · , · , u n ) to g α ( · , · , u ) − g β ( · , · , u ) in L 1 ( S × Ω). F or the second in tegral, Theorem 3.2 sho ws that u n i ( s, · ) → u i ( s, · ) in H 1 0 (Ω) for each i and almost surely in s ∈ S . As in the pro of of Theorem 2.5 , using the a priori b ound ( 9 ) and the dominated con vergence theorem, we obtain the strong conv ergence u n → u in L p (Ω; H 1 0 (Ω; R m )), i.e. ∇ u n i → ∇ u i in L p (Ω; L 2 (Ω; R n )). T o obtain an analogous conclusion for the adjoin t functions p n , w e first show the strong con vergence ∂ g γ ∂ u i ( · , · , u n ) → ∂ g γ ∂ u i ( · , · , u ) in L p ′ ( S ; L q ′ ∗ (Ω)) for γ ∈ { α, β } . (18) This can b e prov ed in a manner analogous to the pro of of Theorem 2.5 . The only difference o ccurs in the step following ( 12 ), where w e ma y tak e q ∗ instead of q ∈ [1 , q ∗ ), since the compactness of the em b edding H 1 0 (Ω)  → L q (Ω) is not required. As a result, u n → u strongly in L p ( S ; L q ∗ (Ω; R m )), and passing to a subsequence, | u n | ≤ V almost surely-almost everywhere on S × Ω for some V ∈ L p ( S ; L q ∗ (Ω)). Now, along a further subsequence, we ha ve p oin twise a.s.-a.e. con v ergence of the conv ergence that app ears in ( 18 ), and by the gro wth conditions ( 15 ),     ∂ g γ ∂ u i ( s, x , u n ( s, x ))     ≤ e ϕ γ ( s, x ) + e ψ γ ( s, x ) V ( s, x ) q 1 − 1 . 13 W e ha ve V q 1 − 1 ∈ L p q 1 − 1 ( S ; L q ∗ q 1 − 1 (Ω)), so by H¨ older’s inequality e ψ γ V q 1 − 1 ∈ L p ′ ( S ; L q ′ ∗ (Ω)) pro vided that 1 q 4 + q 1 − 1 p = 1 p ′ = p − 1 p , 1 q 5 + q 1 − 1 q ∗ = 1 q ′ ∗ = q ∗ − 1 q ∗ whic h hold by assumption ( 15 ). Applying the dominated conv ergence theo- rem then yields ( 18 ). Then, the righ t-hand side in the adjoin t equation ( − div( A n ∇ p n i ) = θ n ∂ g α ∂ u i ( · , · , u n ) + (1 − θ n ) ∂ g β ∂ u i ( · , · , u n ) p n i ∈ L p ′ ( S ; H 1 0 (Ω)) i = 1 , . . . , m, b elongs to L p ′ ( S ; L q ′ ∗ (Ω)), which is contin uously em b edded in L p ′ ( S ; H − 1 (Ω)), yielding the unique solution p n i ∈ L p ′ ( S ; H 1 0 (Ω)). Moreov er, from ( 18 ) and Theorem 3.2 , we hav e conv ergence p n ( s, · ) → p ( s, · ) in H 1 0 (Ω) almost surely in s . As sho wn for the sequence ( u n ) in the proof of Theorem 2.5 , this implies the strong con vergence p n → p in L p ′ (Ω; H 1 0 (Ω; R m )). Finally , we conclude that ∇ u n i · ∇ p n i → ∇ u i · ∇ p i strongly in L 1 ( S × Ω), whic h completes the pro of. T o write do wn the optimalit y conditions for ( P ), we need some further assumptions on  . One of those requires the notion of r e gularity in the sense of Clark e [ 9 , Definition 2.3.4, p. 39]. Indeed, the function ρ is said to b e r e gular at θ if the (usual one-sided) directional deriv ativ e at θ exists for all directions v and it equals the gener alize d dir e ctional derivative ρ ◦ ( θ ; v ) = lim sup η → θ,t ↓ 0 ( ρ ( η + tv ) − ρ ( η )) /t. W e mak e the following assumption. Assumption 3.4. (i) L et  : L ∞ (Ω; [0 , 1]) → R b e Hadamar d dir e ctional ly differ entiable at θ ∗ . (ii) L et  b e lo c al ly Lipschitz and r e gular at θ ∗ . Let us try to write down a first optimality condition. W e argue in a similar fashion to [ 3 ]. W e start b y recalling the definition of the tangen t (Bouligand) cone of A at the p oin t ( θ , A ): T A ( θ , A ) := { ( h, d ) ∈ L ∞ (Ω; [0 , 1]) × L ∞ (Ω; Sym d ) | ∃ s k ↘ 0 , ∃ ( h k , d k ) → ( h, d ) , ( θ , A ) + s k ( h k , d k ) ∈ A for eac h k } . Lemma 3.5. L et The or ems 2.3 , 2.4 and 3.1 and The or em 3.4 (i) hold. If ( θ ∗ , A ∗ ) is an optimal r elaxe d solution of ( P ) , then, for every ( δ θ, δ A ) ∈ T A ( θ ∗ , A ∗ ) we have sup ν ∈ ∂ R ( J ( θ ∗ , A ∗ )) E [ ν J ′ (( θ ∗ , A ∗ ))( δ θ , δ A )] +  ′ ( θ ∗ )( δ θ ) ≥ 0 . 14 Pr o of. Firstly , since R is finite, low er semicon tinuous and con vex, R is con- tin uous [ 8 , Proposition 2.111], and by [ 8 , Proposition 2.126 (v)], we obtain that in fact R is Hadamard directionally differentiable with R ′ [ z ]( h ) = sup ν ∈ ∂ R ( z ) E [ ν h ] . No w we can use the c hain rule [ 8 , Prop osition 2.47]: since J : A → L 1 (Ω) is C 1 b y Theorem 3.3 and hence also Hadamard directionally differentiable, and R is Hadamard directionally differentiable, w e find that the comp osition R ◦ J is Hadamard directionally differentiable to o with ( R ◦ J ) ′ ( θ ∗ , A ∗ )( δ θ , δ A ) = R ′ [ J ( θ ∗ , A ∗ )] J ′ (( θ ∗ , A ∗ ))( δ θ , δ A ) , for ev ery ( δ θ , δ A ) ∈ T A ( θ ∗ , A ∗ ). The set A ⊂ L ∞ (Ω; [0 , 1]) × L ∞ (Ω; Sym d ) is in general non-con vex, hence the closure of the radial cone ma y not b e the tangen t cone, so the usual density argumen t to obtain the claim of this result has to be adapted. No w, if ( δ θ , δ A ) ∈ T A ( θ ∗ , A ∗ ), by definition, we ha ve the existence of { ( h k , d k ) } suc h that ( h k , d k ) → ( δ θ , δ A ) and a n ull sequence { s k } with R [ J ( θ ∗ + s k h k , A ∗ + s k d k )] +  ( θ ∗ + s k h k ) − R [ J ( θ ∗ , A ∗ )] −  ( θ ∗ ) ≥ 0 , since ( θ ∗ , A ∗ ) is optimal. Dividing the ab o ve expression b y s k and taking the limit k → ∞ and using the fact that R ◦ J and  are Hadamard directionally differen tiable, we get R ′ [ J ( θ ∗ , A ∗ )]( J ′ (( θ ∗ , A ∗ ))( δ θ , δ A )) +  ′ ( θ ∗ )( δ θ ) ≥ 0 . Using the characterisation of R ′ in terms of the exp ectation stated ab o ve, w e obtain the result. As done in [ 18 ], w e wish to reform ulate the ab o v e condition since the presence of the suprem um renders it inconv enien t and mak es it less tractable. By means of sub differen tial calculus, this reformulation is p ossible as w e see in the next result. This result enables the adaptation of the optimality criteria metho d which is well established in the deterministic setting to the risk- a verse setting as well. Belo w, ∂  refers to the generalized gradient in the sense of Clark e, see [ 9 , p. 27]. Theorem 3.6. L et ( θ ∗ , A ∗ ) b e an optimal r elaxe d solution of ( P ) . L et The o- r ems 2.3 , 2.4 , 3.1 and 3.4 hold. Then, ther e exists π ∗ ∈ ∂ R ( J ( θ ∗ , A ∗ )) such that for every ( δ θ , δ A ) ∈ T A ( θ ∗ , A ∗ ) , 0 ≤  ′ ( θ ∗ )( δ θ ) + Z S Z Ω π ∗ δ θ  g α ( s, x , u ∗ ) − g β ( s, x , u ∗ )  d x d P ( s ) − Z S Z Ω π ∗ m X i =1 δ A ∇ u ∗ i · ∇ p ∗ i d x d P ( s ) . (19) 15 Pr o of. The pro of is similar to that of [ 3 , Lemma 3.11] with some adaptations to deal with the non-conv exity . T o b egin, observ e that R is lo cally Lipschitz near every p oin t of L 1 (Ω) and that J : A → L 1 (Ω), b eing C 1 , is also strictly differen tiable (in the sense of Clarke) [ 9 , Corollary , p. 32]. By [ 9 , Prop osition 2.2.1], J is Lipschitz near ( θ ∗ , A ∗ ). Therefore, the comp osition R ◦ J is lo cally Lipschitz near ( θ ∗ , A ∗ ). No w,  is assumed lo cally Lipsc hitz on L ∞ (Ω; [0 , 1]). In com bination with the ab o ve, w e thus ha ve ha v e that the sum R ◦ J +  is also lo cally Lipsc hitz. In voking the corollary on page 52 of [ 9 ], w e get 0 ∈ ∂ ( R ◦ J +  )( θ ∗ , A ∗ ) + N A ( θ ∗ , A ∗ ) where N A ( u ) = T A ( u ) ◦ = { η ∈ (L ∞ (Ω; [0 , 1]) × L ∞ (Ω; Sym d )) ∗ : ⟨ η , v ⟩ ≤ 0 , v ∈ T A ( u ) } is the normal cone of A at the p oin t u . Applying [ 9 , Prop o- sition 2.3.3], the first term on the right-hand side ab ov e can b e expressed as ∂ ( R ◦ J +  )( θ ∗ , A ∗ ) ⊂ ∂ ( R ◦ J )( θ ∗ , A ∗ ) + ∂  ( θ ∗ ) . Since R is locally Lipsc hitz on L 1 (Ω), w e are able to utilize the sub differential c hain rule [ 9 , Theorem 2.3.10 and Remark 2.3.11] to find ∂ ( R ◦ J )( θ ∗ , A ∗ ) = [ J ′ ( θ ∗ , A ∗ )] ∗ ∂ R ( J ( θ ∗ , A ∗ )) (equalit y holds since R is regular [ 9 , Proposition 2.3.6 (b)] by virtue of its con vexit y). Gathering ev erything together, w e ha ve 0 ∈ [ J ′ ( θ ∗ , A ∗ )] ∗ ∂ R ( J ( θ ∗ , A ∗ )) + ∂  ( θ ∗ ) + N A ( θ ∗ , A ∗ ) . No w, as in the pro of of [ 3 , Lemma 3.11], w e follo w [ 18 ]. F rom the ab o v e inclusion, w e obtain the existence of π ∗ ∈ ∂ R ( J ( θ ∗ , A ∗ )) and η ∈ ∂  ( θ ∗ ) suc h that − [ J ′ ( θ ∗ , A ∗ )] ∗ π ∗ − η ∈ N A ( θ ∗ , A ∗ ) . This reads ⟨− [ J ′ ( θ ∗ , A ∗ )] ∗ π ∗ − η , v ⟩ ≤ 0 , v ∈ T A ( θ ∗ , A ∗ ) . No w, since  is regular, b y definition of the generalized gradien t, we hav e that ∂  ( θ ∗ ) = { ξ ∈ (L ∞ (Ω; [0 , 1])) ∗ :  ′ ( θ ∗ )( v ) ≥ ⟨ ξ , v ⟩ , v ∈ L ∞ (Ω; [0 , 1]) } . Th us the ab o v e can b e written as ⟨− [ J ′ ( θ ∗ , A ∗ )] ∗ π ∗ −  ′ ( θ ∗ ) , v ⟩ ≤ 0 , v ∈ T A ( θ ∗ , A ∗ ) . 16 Then E [ J ′ ( θ ∗ , A ∗ ) v π ∗ ] +  ′ ( θ ∗ )( δ θ ) ≥ 0 , v = ( δ θ , δ A ) ∈ T A ( θ ∗ , A ∗ ) follo ws up on realizing ⟨J ′ ( θ ∗ , A ∗ ) ∗ π ∗ , v ⟩ = ⟨ π ∗ , J ′ ( θ ∗ , A ∗ ) v ⟩ as the L 1 ( S )- L 1 ( S ) ∗ pairing. The claim of the theorem is a direct consequence of substi- tuting the deriv ativ e form ula ( 16 ). T o further analyze the optimalit y condition of Theorem 3.6 , let us sp e- cialize to the case where the regularization term is of the form  ( θ ) = Z Ω h ( θ ) d x , (20) where h : [0 , 1] → R is a con tinuously differentiable function. In classical applications, a choice suc h as  ( θ ) = λ R Ω θ d x p enalizes the amoun t of the first phase. F or greater generalit y , w e consider instead the form stated ab ov e. W e can analyze the result of Theorem 3.6 using the same approac h as in the deterministic setting [ 22 , 1 ]. The dep endence on the v ariable s in the optimalit y conditions can be addressed in a manner reminiscent of the treatmen t of the time v ariable in evolutionary problems [ 31 ]. W e b egin b y fixing θ and considering v ariations only in A . Since K ( ϑ ) is conv ex for every ϑ ∈ [0 , 1], then for every A ∈ K ( θ ∗ ) w e hav e (0 , A − A ∗ ) ∈ T A ( θ ∗ , A ∗ ). Therefore, w e obtain − Z S Z Ω π ∗ m X i =1 ( A − A ∗ ) ∇ u ∗ i · ∇ p ∗ i d x d P ( s ) ≥ 0 . By in terchanging the order of in tegration, this ma y b e rewritten as − Z Ω Z S π ∗ m X i =1 ( A − A ∗ ) ∇ u ∗ i · ∇ p ∗ i d P ( s ) d x ≥ 0 . W e define M ∗ := Sym Z S π ∗ m X i =1 ∇ u ∗ i ⊗ ∇ p ∗ i d P ( s ) , (21) where Sym denotes the symmetric part of a matrix, and ⊗ the tensor pro duct of t wo vectors in R d . With this notation, the ab o ve condition can b e rewritten as Z Ω A ∗ : M ∗ d x ≥ Z Ω A : M ∗ d x , 17 where : denotes the Euclidean inner pro duct on M d ( R ). Since A ∈ K ( θ ∗ ) is arbitrary , we may choose A to equal A ∗ a.e. except in a small neighborho o d of an arbitrary point. This yields the following p oin twise estimate: almost ev erywhere on Ω, A ∗ : M ∗ ≥ A : M ∗ for an y A ∈ K ( θ ∗ ). T o summa- rize, the necessary optimality condition can b e expressed as follo ws: almost ev erywhere on Ω, the optimal A ∗ satisfies A ∗ : M ∗ = F ( θ ∗ , M ∗ ) , where F : [0 , 1] × M d ( R ) → R is giv en by F ( ϑ, M ) = max A ∈ K ( ϑ ) A : M . W e can now pro ceed as in [ 1 , Theorem 3.2.14] and consider v ariations in θ . Since the function F is con tin uously differen tiable with resp ect to θ , w e can consider a smo oth path θ ( t ) ∈ L ∞ (Ω; [0 , 1]) passing through θ ∗ with deriv ativ e δθ , and c ho ose A ( t ) ∈ K ( θ ( t )) such that A ( t ) : M ∗ = F ( θ ( t ) , M ∗ ), almost everywhere on Ω. Consequently , the deriv ativ e of A with resp ect to t , ev aluated at t = 0 where θ (0) = θ ∗ , satisfies δ A : M ∗ = ∂ F ∂ θ ( θ ∗ , M ∗ ) δ θ . The partial deriv ativ e of F with resp ect to θ can b e expressed explicitly in the t wo- and three-dimensional case [ 4 , 30 ]. Finally , w e come to the following result. Theorem 3.7. L et  b e of the form ( 20 ) and let the assumptions of The o- r em 3.6 hold. If ( θ ∗ , A ∗ ) is an optimal r elaxe d solution of ( P ) and u ∗ and p ∗ the c orr esp onding state and adjoint state functions r esp e ctively, then, ther e exists π ∗ ∈ ∂ R ( J ( θ ∗ , A ∗ )) such that for Q : Ω → R define d by Q = Z S π ∗  g α ( s, · , u ∗ ) − g β ( s, · , u ∗ )  d P ( s ) − ∂ F ∂ θ ( θ ∗ , M ∗ ) + h ′ ( θ ) , the fol lowing ne c essary optimality c ondition hold: for almost every x ∈ Ω Q ( x ) > 0 = ⇒ θ ∗ ( x ) = 0 Q ( x ) < 0 = ⇒ θ ∗ ( x ) = 1 . 4 Numerics via the optimalit y criteria metho d In the section, we will study a particular instance of ( P ) and solve it numer- ically . W e will tak e adv an tage of the first-order conditions we deriv ed in the previous section and use the optimal criteria method (see e.g. [ 5 , Chapter 1] 18 for more details) in order to compute a solution of this problem. The metho d in our general problem setting is given in Algorithm 4.1 . Algorithm 4.1 Optimality Criteria Metho d Initialization. Cho ose an admissible initial design ( θ 0 , A 0 ). F or k = 0 , 1 , 2 , . . . , rep eat the following steps: Step 1 (State equation). Compute the state u k = ( u k 1 , . . . , u k m ) = G ( A k ) . Step 2 (Adjoin t equation). Compute the adjoint v ariable p k = ( p k 1 , . . . , p k m ) , solving ( 17 ) with A k , θ k , and u k in place of θ ∗ , A ∗ , and u ∗ . Step 3 (Sensitivit y matrix). Compute M k = Sym Z S π k m X i =1 ∇ u k i ⊗ ∇ p k i d P ( s ) , where π k is a subgradien t of R at J ( θ k , A k ). Step 4 (Up date of θ ). F or each x ∈ Ω, define θ k +1 ( x ) as a zero of θ 7→ Z S π k  g α ( s, x , u k ) − g β ( s, x , u k )  d P ( s ) − ∂ F ∂ θ ( θ , M k ( x )) + h ′ ( θ ) . (22) If a zero do es not exist, set θ k +1 ( x ) = 0 when the function is p ositive, and θ k +1 ( x ) = 1 when it is negativ e. Step 5 (Up date of A). F or each x ∈ Ω, choose A k +1 ( x ) ∈ arg max A F  θ k +1 ( x ) , M k ( x )  . W e still need to sp ecify for example ho w to numerically appro ximate the in tegrals o ver the probabilit y space S that app ear in the algorithm; this will b e addressed in the next section for our sp ecific example. 4.1 An example in volving the conditional v alue-at-risk Let us now sp ecify the exact problem w e will simulate. F or the risk measure R , w e adopt the conditional (a v erage) v alue-at-risk CV aR. T o define this, 19 let us recall some definitions from [ 24 , § 6.2.3–6.2.4] for the conv enience of the reader. Firstly , the value-at-risk (V aR) at lev el γ ∈ (0 , 1) for a random v ariable representing losses, also known as the γ -quan tile, is defined as V aR γ ( X ) = inf { t : P ( X > t ) ≤ γ } = inf { t : P ( X ≤ t ) ≥ 1 − γ } . Throughout this text, we assume that the random v ariable X has no proba- bilit y atoms, i.e. its distribution is absolutely contin uous with resp ect to the Leb esgue measure. The c onditional value-at-risk (CV aR) at level γ is the exp ected loss given that the loss exceeds the v alue-at-risk V aR γ at lev el γ : CV aR γ ( X ) = E [ X | X ≥ V aR γ ( X )] . (23) In tuitively , CV aR γ can b e understoo d as the exp ected loss in the w orst γ - fraction or tail of outcomes. T o keep the notation simple, w e consider the case of a single state problem (i.e., m = 1)  − div( A ∇ u ) = f ( s, x ) u ∈ H 1 0 (Ω) , (24) where Ω = [0 , 1] 2 is the square. In the context of heat conduction, u denotes the temp erature field, f represents the densit y of the external heat source, and the ob jective is to maximize the total heat energy R Ω f u d x . Therefore, in our risk-a verse formulation, we seek to minimize large deviations of J ( θ , A ) = − Z Ω u ( s, x ) f ( s, x ) d x . The function  represents the constrain t on the total amount of the first phase, giv en by  ( θ ) = λ Z Ω θ , where λ is the asso ciated Lagrange m ultiplier. Note that the sub differen tial of R = CV aR γ can b e characterized as follo ws. F or a giv en X ∈ L 1 ( S ), the subdifferential ∂ CV aR γ ( X ) is a singleton con taining the unique element π ∈ L ∞ ( S ) defined b y [ 24 ] π ( s ) = ( 0 , X ( s ) < V aR γ ( X ) 1 γ , X ( s ) > V aR γ ( X ) (b earing in mind that P ( X = V aR γ ( X )) = 0). T o facilitate the ev aluation of V aR γ ( J ( θ , A )), w e employ the Karh unen–Lo ` ev e sp ectral decomp osition, whic h separates the spatial and sto c hastic comp onen ts: f ( s, x ) = f 0 ( x ) + ∞ X j =1 p λ j f j ( x ) ξ j ( s ) , 20 where { ξ i } i ∈ N are mutually uncorrelated random v ariables with mean zero and unit v ariance. W e fo cus on the truncated sum f ( s, x ) = f 0 ( x ) + N X j =1 p λ j f j ( x ) ξ j ( s ) . If u j denotes the solution of the deterministic state equation ( 6 ) with the righ t-hand side f j , for j = 0 , . . . , N , then the p erturb ed state function is simply giv en by u ( s, x ) = u 0 ( x ) + N X j =1 p λ j u j ( x ) ξ j ( s ) , whic h leads to J ( θ , A )( s ) = c 0 + N X j =1 c j ξ j ( s ) + N X i =1 N X j = i c ij ξ i ( s ) ξ j ( s ) where c 0 = − Z Ω u 0 ( x ) f 0 ( x ) d x , c j = − p λ j Z Ω u 0 ( x ) f j ( x ) + u j ( x ) f 0 ( x ) d x , j = 1 , . . . , N , c ij = − p λ i λ j Z Ω u i ( x ) f j ( x ) + u j ( x ) f i ( x ) d x , i, j = 1 , . . . , N , j ≥ i . F or the optimalit y criteria metho d, let us denote V aR γ ( J ( θ k , A k )) simply b y δ k , and in tro duce the set ∆ k = { s ∈ S : J ( θ k , A k )( s ) ≥ δ k } , so that the subgradient π k is giv en by 1 γ χ ∆ k . In the problem under consideration, w e hav e g α = g β , which implies that the first in tegral in ( 22 ) v anishes. F urthermore, the problem is self-adjoin t; more precisely p k = − u k , so that M k = −∇ u k 0 ⊗ ∇ u k 0 − 2 N X j =1 d k j Sym d ∇ u k 0 ⊗ ∇ u k j − 2 N X i =1 N X j = i d k ij Sym d ∇ u k i ⊗ ∇ u k j where d k j = 1 γ p λ j Z ∆ k ξ j d P ( s ) , j = 1 , . . . , N , d k ij = 1 γ p λ i λ j Z ∆ k ξ i ξ j d P ( s ) , i, j = 1 , . . . , N , j ≥ i . 21 W e employ an expansion based on an exp onen tial co v ariance function. W e b egin with the zero-mean random field g ( s, x ) = f ( s, x ) − f 0 ( x ), since E [ f ( · , x )] = f 0 ( x ). The cov ariance function is given by C ( x , x ′ ) = E [ g ( · , x ) g ( · , x ′ )] , and in this example we sp ecify C to b e exp onen tial cov ariance function: C ( x , x ′ ) = σ 2 exp − d X k =1 | x k − x ′ k | η k ! , where σ 2 and η k are the v ariance and the correlation lengths. F or this c hoice of cov ariance, the eigenpairs ( λ j , f j ) of the associated in tegral operator C : L 2 (Ω) → L 2 (Ω) defined b y C ψ ( x ) = Z Ω C ( x , x ′ ) ψ ( x ′ ) d x ′ can be computed explicitly [ 35 ]. Moreov er, the random co efficien ts ξ j in the resulting expansion are standard normally distributed. The n umerical example considers the tw o-dimensional setting Ω = [0 , 1] 2 (as mentioned b efore), with a constan t mean v alue f 0 = 1 on Ω. W e choose N = 5, σ = 2, lev el γ = 0 . 9, and examine sev eral v alues of the correlation length η := η 1 = η 2 . The algorithm is implemented in Freefem++ , with P1 elemen ts. A t each iteration of the optimality criteria method, the γ -quantile V aR γ ( J ( θ k , A k )) and the in tegrals app earing in the definition of M k are appro ximated using a Monte Carlo pro cedure with a sample size of 10,000. As is typical of the optimality criteria metho d, only a few iterations are required to obtain a go o d approximation of the solution, and the result is largely insensitiv e to the c hoice of the initial design. Figure 1 presents the designs obtained after 10 iterations for three v alues of the parameter η , along- side the design corresp onding to the unp erturb ed right-hand side f 0 = 1 for comparison. 22 Figure 1: Numerical solutions for η = 0 . 5, 0 . 1, and 0 . 05 (display ed left to righ t, top to b ottom) for R = CV aR γ , along with the solution corresp onding to the unperturb ed right-hand side. W e see that as η shrinks (i.e. as the influence of randomness reduces), the optimal shap es resemble more and more the deterministic shap e seen in the b ottom right corner of Figure 1 . F or purp oses of comparison, we also consider the risk-neutral case in whic h the risk measure coincides with the exp ectation R = E . The corre- sp onding designs are sho wn in Figure 2 , again for v arious v alues of η and the unp erturbed righ t-hand side in the rightmost panel. It can b e observed that the obtained solutions are close to those of the corresp onding non-p erturb ed problem. This b eha vior can b e explained by the fact that the random fluc- tuations ha ve zero mean, and th us tend to cancel out in exp ectation. In contrast, for the risk-av erse form ulation, extreme realizations dominate the ob jective, and the optimizer explicitly accoun ts for adv erse ev ents. As a result, the solution b ecomes more conserv ative, prioritizing robustness o ver a verage performance. This approach reduces the lik eliho o d of po or outcomes under unfav orable realizations, even at the exp ense of sligh tly sub optimal p erformance under typical conditions. 23 Figure 2: (Risk-neutral case) Numerical solutions for η = 0 . 5, 0 . 05 (displa yed left to righ t) for R = E , along with the solution corresponding to the unper- turb ed right-hand side. 5 Conclusion In this pap er, we studied a risk-a verse stochastic optimization problem in- v olving optimal designs. W e prov ed existence and deriv ed stationarity con- ditions, and used those to set up a n umerical scheme and n umerically solved an example. Our work can be seen as a first step in taking into account risk measures to deal with uncertain ties in optimal design problems in conduc- tivit y . F urther work could inv olve a more sophisticated analysis sp ecialized to certain problem classes in applications and the resulting n umerics and/or the dev elopment of tailored solution algorithms. Ac kno wledgemen ts The research of P . Kun ˇ stek and M. V rdoljak has b een supp orted in part b y Croatian Science F oundation under the pro ject IP-2022-10-7261. The calcu- lations were p erformed at the Lab oratory for Adv anced Computing (F aculty of Science, Universit y of Zagreb). A. Alphonse thanks Ren ´ e Henrion for very helpful discussions. References [1] G. Allaire. Shap e optimization by the homo genization metho d , v olume 146 of Applie d Mathematic al Scienc es . Springer-V erlag, New Y ork, 2002. [2] G. Allaire and C. Dapogny . A deterministic approximation metho d 24 in shap e optimization under random uncertain ties. SMAI J. Comput. Math. , 1:83–143, 2015. [3] A. Alphonse, C. Geiersbach, M. Hinterm¨ uller, and T. M. Surowiec. Risk- a verse optimal con trol of random elliptic v ariational inequalities. ESAIM Contr ol Optim. Calc. V ar. , 31, 2025. P ap er No. 71, 38. [4] N. Antoni ´ c and M. V rdoljak. Sequential laminates in m ultiple state optimal design problems. Math. Pr obl. Eng. , 2006:Article ID 68695, 2006. [5] M. P . Bendsøe and O. Sigmund. T op olo gy optimization . Springer-V erlag, Berlin, 2003. [6] M. P . Bendsøe and N. Kikuc hi. Generating optimal top ologies in struc- tural design using a homogenization metho d. Computer Metho ds in Applie d Me chanics and Engine ering , 71(2):197–224, 1988. [7] A. Benedek and R. Panzone. The space L p , with mixed norm. Duke Mathematic al Journal , 28(3):301 – 324, 1961. [8] J. F. Bonnans and A. Shapiro. Perturb ation Analysis of Optimization Pr oblems . Springer Science & Business Media, 2013. [9] F. H. Clarke. Optimization and Nonsmo oth A nalysis , v olume 5 of Clas- sics in Applie d Mathematics . So ciet y for Industrial and Applied Math- ematics (SIAM), Philadelphia, P A, second edition, 1990. [10] S. Conti, H. Held, M. P ach, M. Rumpf, and R. Sc h ultz. Risk av erse shap e optimization. SIAM J. Contr ol Optim. , 49(3):927–947, 2011. [11] M. Dambrine, H. Harbrec h t, and B. Puig. Computing quan tities of in terest for random domains with second order shap e sensitivit y analysis. ESAIM Math. Mo del. Numer. A nal. , 49(5):1285–1302, 2015. [12] H. Ga jewski, K. Gr¨ oger, and K. Zac harias. Nichtline ar e Op er ator gle- ichungen und Op er ator differ entialgleichungen . Ak ademie-V erlag, Berlin, 1974. [13] A. Hantoute, R. Henrion, and P . P´ erez-Aros. Sub differen tial c harac- terization of probabilit y functions under Gaussian distribution. Math. Pr o gr am. , 174(1-2):167–194, 2019. [14] M. Heink ensc hloss and D. P . Kouri. Optimization problems go v erned by systems of PDEs with uncertainties. A cta Numer. , 34:491–577, 2025. 25 [15] T. Hyt¨ onen, J. v an Neerv en, M. V eraar, and L. W eis. Analysis in Banach Sp ac es: V olume I: Martingales and Littlewo o d-Paley The ory . Springer, 2016. [16] D. P . Kouri and T. M. Suro wiec. Existence and Optimalit y Conditions for Risk-Av erse PDE-Constrained Optimization. SIAM/ASA Journal on Unc ertainty Quantific ation , 6(2):787–815, 2018. [17] D. P . Kouri and T. M. Surowiec. Risk-a v erse optimal control of semilin- ear elliptic PDEs. ESAIM Contr ol Optim. Calc. V ar. , 26, 2020. P ap er No. 53, 19. [18] D. P . Kouri and T. M. Suro wiec. Corrigendum: “Existence and optimalit y conditions for risk-av erse PDE-constrained optimization”. SIAM/ASA J. Unc ertain. Quantif. , 10(3):1321–1322, 2022. [19] A. J. Kurdila and M. Zabarankin. Convex functional analysis . Systems & Control: F oundations & Applications. Birkh¨ auser V erlag, Basel, 2005. [20] K. Lurie and A. Cherk aev. Exact estimates of conductivity of comp os- ites formed b y tw o isotropically conducting madia, tak en in prescrib ed prop ortion. Pr o c. R oyal So c. Edinbur gh , 99A:71–87, 1984. [21] J. Martinez-F rutos and F. Periago. Optimal Contr ol of PDEs Under Unc ertainty: An intr o duction with Applic ation to Optimal Shap e Design of Structur es . Springer, 2018. [22] F. Murat and L. T artar. Calcul des v ariations et homog ´ en ´ eisation. In Homo genization metho ds: the ory and applic ations in physics (Br´ eau- sans-Napp e, 1983) , v olume 57 of Col le ct. Dir. ´ Etudes R e ch. ´ Ele c. F r anc e , pages 319–369. Eyrolles, Paris, 1985. [23] F. Murat and L. T artar. H -conv ergence. In T opics in the mathematic al mo del ling of c omp osite materials , volume 31 of Pr o gr. Nonline ar Differ- ential Equations Appl. , pages 21–43. Birkh¨ auser Boston, Boston, MA, 1997. [24] A. Shapiro, D. Dentc hev a, and A. Ruszczy ´ nski. L e ctur es on Sto chas- tic Pr o gr amming: Mo deling and The ory, Se c ond Edition . So ciet y for Industrial and Applied Mathematics, 2014. [25] O. Sigmund and K. Maute. T op ology optimization approaches: A comparativ e review. Structur al and Multidisciplinary Optimization , 48(6):1031–1055, 2013. 26 [26] L. T artar. An introduction to the homogenization metho d in optimal design. In Optimal shap e design (Tr´ oia, 1998) , v olume 1740 of L e ctur e Notes in Math. , pages 47–156. Springer, Berlin, 2000. [27] W. v an Ac kooij and R. Henrion. Gradient form ulae for nonlinear proba- bilistic constraints with Gaussian and Gaussian-lik e distributions. SIAM J. Optim. , 24(4):1864–1889, 2014. [28] W. v an Ack o oij, R. Henrion, and H. Zidani. P on try agin’s principle for some probabilistic control problems. Appl. Math. Optim. , 90(1):P ap er No. 5, 36, 2024. [29] W. v an Ack o oij and P . P´ erez-Aros. Gradien t form ulae for nonlinear probabilistic constraints with non-conv ex quadratic forms. J. Optim. The ory Appl. , 185(1):239–269, 2020. [30] M. V rdoljak. On Hashin-Sh trikman b ounds for mixtures of t w o isotropic materials. Nonline ar Anal. R e al World Appl. , 11(6):4597–4606, 2010. [31] M. V rdoljak. Optimalit y criteria metho d for optimal design in hyperb olic problems. Math. Commun. , 15:555–572, 2010. [32] M. V rdoljak. Classical optimal design in t wo-phase conductivity prob- lems. SIAM J. Contr ol Optim. , 54(4):2020–2035, 2016. [33] M. V rdoljak. Multiple state optimal design problems with random p er- turbation. Ele ctr on. J. Differ ential Equations , 2018. P ap er No. 59, 10. [34] K. Y osida. F unctional analysis . Classics in Mathematics. Springer- V erlag, Berlin, 1995. Reprint of the sixth (1980) edition. [35] D. Zhang and Z. Lu. An efficient, higher-order p erturbation approach for flo w in randomly heterogeneous p orous media via Karh unen–Lo` ev e de- comp osition. Journal of Computational Physics , 192:773–794, 01 2004. 27