\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2018 (2018), No. 59, pp. 1--10.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2018 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2018/59\hfil 
Multiple state optimal design problems]
{Multiple state optimal design problems with random perturbation}

\author[M. Vrdoljak \hfil EJDE-2018/59\hfilneg]
{Marko Vrdoljak}

\address{Marko Vrdoljak \newline
Department of Mathematics,
Faculty of Science,
University of Zagreb, Croatia}
\email{marko@math.hr}

\thanks{Submitted April 28, 2017. Published March 2, 2018.}
\subjclass[2010]{49K35, 49K20, 49J20, 80M40}
\keywords{Stationary diffusion; optimal design; homogenization;
 \hfill\break\indent random perturbation; optimality conditions}

\begin{abstract}
 A multiple state optimal design problem with  presence of uncertainty on
 the right-hand side is considered, in the context of  stationary diffusion
 with two isotropic  phases.
 A similar problem with one state equation has already been considered by
 Buttazzo and Maestre (2011). We shall address the question of relaxation
 by the homogenization method and necessary conditions of optimality.
 The case of discrete probability space  leads to another multiple state problem
 (possibly with an infinite number of states), which could be treated by similar
 techniques to those presented in Allaire (2002) and Vrdoljak (2010).
 The relaxation can be expressed in a simpler form for problems with spherical
 symmetry in the case of minimization (or maximization) of averaged
 energy, and we present an example which can be solved explicitly. 
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

For a measurable matrix function $\mathbf{A}$ which is bounded and uniformly elliptic
 (almost everywhere on an open and bounded set $\Omega\subseteq\mathbb{R}^d$)  and 
$f\in H^{-1}(\Omega)$ there exists a unique solution of the
boundary value problem for stationary diffusion equation
\begin{equation}\label{sdif}
\begin{gathered}
-\operatorname{div}(\mathbf{A}\nabla u)=f\\
u\in {\rm H}^1_0(\Omega)\,.
\end{gathered}
\end{equation}

The homogenization theory allows one to introduce a topology on the appropriate 
set of coefficients such that the mapping $f\mapsto u$ is continuous, in a 
reasonable pair of topologies. Historically, these topologies were first 
introduced (with full mathematical rigour) by Spagnolo \cite{Ssemel}
through the concept of $G$-convergence.
The notion of $H$-convergence was also originally introduced for the stationary 
diffusion equation \cite{MT1}
(it is also known under the name strong G-convergence \cite{ZKON}).

More precisely, Murat and Tartar \cite{MT1} introduced the set of admissible
conductivity matrix functions
\begin{align*}
\mathcal{M}(\alpha,\beta;\Omega)=\Big\{&\mathbf{A}\in L^\infty(\Omega; M_d(\mathbb{R})): 
 \mathbf{A}(\mathbf{x})\boldsymbol{\xi}\cdot\boldsymbol{\xi}\geq\alpha |\boldsymbol{\xi}  |^2\,,	\\
&\mathbf{A}(\mathbf{x})\boldsymbol{\xi}\cdot\boldsymbol{\xi}\geq {\frac1\beta}|\mathbf{A}(\mathbf{x})\boldsymbol{\xi}|^2\Big\}\,.
\end{align*}
A sequence of matrix functions $(\mathbf{A}_n)$ in $\mathcal{M}(\alpha,\beta;\Omega)$ is said to $H$-converge 
to $\mathbf{A}\in\mathcal{M}(\alpha',\beta';\Omega)$ if for any $f\in H^{-1}(\Omega)$ the sequence $(u_n)$ of
 solutions of \eqref{sdif}, with $\mathbf{A}_n$ instead of $\mathbf{A}$, satisfies
\begin{gather}\label{H}
u_n \rightharpoonup u \quad \text{in }  H_0^1(\Omega)\,,\\
\mathbf{A}_n\nabla u_n \rightharpoonup \mathbf{A} \nabla u \quad \text{in }  L^2(\Omega;\mathbb{R}^d)\,,
\end{gather}
where $u$ is given by \eqref{sdif}.
If we additionally assume that conductivity matrices are symmetric,
 the above definition coincides with the notion of $G$-convergence,
imposing only the convergence \eqref{H}.
The question whether such $H$-converging sequence exists is simply answered 
by the compactness theorem: the bounds in the
definition of set $\mathcal{M}(\alpha,\beta;\Omega)$ are well chosen in such a way that it is compact 
with respect to $H-$convergence. In other words, in
the definition of $H$-convergence
one can put $\alpha'=\alpha$ and $\beta'=\beta$.

In multiple state optimal design problems, one is trying to find the
best arrangement of given materials, such that the obtained body has
some optimal properties regarding $m$  different regimes.
We shall study the simplest case of two isotropic constituents, with
conductivities $\alpha$ and $\beta$. Therefore, the
conductivity can be written as
$\mathbf{A}=\chi\alpha\mathbf{I}+(1-\chi)\beta\mathbf{I}$,
 where $\chi$ is the characteristic function of the first material.
The optimality of a distribution is measured by an objective function, 
which is usually  an integral functional depending on the distribution of 
materials and the state function, obtained as a solution
of the associated boundary value problem.
If $\chi$ denotes a characteristic function of the first material on a bounded 
open set $\Omega\subseteq \mathbb{R}^d$, this functional can be written 
in the general form
$$
J(\chi)=\int_\Omega \Bigl(\chi(\mathbf{x})\, g_\alpha(\mathbf{x},{\sf u}) 
+ (1-\chi(\mathbf{x})) g_\beta(\mathbf{x},{\sf u})\Bigr)\, d\mathbf{x}\,,
$$
 where ${\sf u}=(u_1,\ldots,u_m)$ is the state function determined by
\begin{equation}\label{state}
\begin{gathered}
-\operatorname{div}(\mathbf{A}\nabla u_i)=f_i\\
u_i\in {\rm H}^1_0(\Omega)
\end{gathered}
\end{equation}
for $i=1,\ldots, m$,
with $\mathbf{A}=\chi\alpha\mathbf{I}+(1-\chi)\beta\mathbf{I}$.  The functions
$f_1,\ldots, f_m$, as well as $g_\alpha$ and $g_\beta$ are supposed to be given.
 The case $m=1$ is studied in \cite{MT,Ta1,Ta2,A} and
the general case in \cite{A,AVmpe,AVfer}.
As it is common in optimal design problems, due to the lack of existence of 
a classical optimal design, one can use the relaxation by the
homogenization method, introducing generalized designs that correspond to 
fine mixtures of original phases.
Under some growth conditions on $g_\alpha$ and $g_\beta$ 
(see \cite[Section 3.1.3]{A}) the relaxed problem reads
\begin{gather*}
J(\theta,\mathbf{A})=\int_\Omega \Bigl(\theta(\mathbf{x})\, g_\alpha(\mathbf{x},{\sf u}(\mathbf{x})) 
+ (1-\theta(\mathbf{x})) g_\beta(\mathbf{x},{\sf u}(\mathbf{x}))\Bigr)\, d\mathbf{x}\to \min\\
\theta\in {\rm L}^\infty(\Omega;[0,1])\,,\;
\mathbf{A}\in\mathcal{K}(\theta) \text{ a.e. on } \Omega\,, \;{\sf u} \text{ solves }\eqref{state}\,.
\end{gather*}

The purpose of this paper is to study  uncertainty perturbation of the 
right-hand side of state equations, contrary to the above deterministic case.
To model such problem we take into consideration a probability space $(S,\mathcal{M},\mu)$.
Suppose that $f_i\in L^1(S,H^{-1}(\Omega))$, and denote
$\overline f_i:=\int_S f_i\,d\mu$.
In other words we consider $s\in S$ to be a parameter in the boundary
 value problem
\begin{equation}\label{ses}
\begin{gathered}
-\operatorname{div}(\mathbf{A}\nabla u_i)=f_i(s,\cdot)\\
u_i\in {\rm H}^1_0(\Omega)\,
\end{gathered}
\end{equation}
for $i=1\ldots m$.

We  consider the following  optimal design problem: 
Given $f_i\in L^1(S,H^{-1}(\Omega))$, $i=1\ldots m$, one seeks for
a  characteristic
function $\chi$ on $\Omega$ that minimizes the following functional
$$
J(\chi)=\int_S\int_\Omega \Bigl(\chi(\mathbf{x})
 g_\alpha(s,\mathbf{x},{\sf u}(s,\mathbf{x})) 
+ (1-\chi(\mathbf{x}))\, g_\beta(s,\mathbf{x},{\sf u}(s,\mathbf{x})\Bigr)\, d\mathbf{x}\,d\mu\,,
$$
where $u_i$ is determined by \eqref{ses} with $\mathbf{A}=\chi\alpha\mathbf{I}+(1-\chi)\beta\mathbf{I}$, 
and ${\sf u}=(u_1,\ldots,u_m)$. Moreover, we assume that the quantity
of the first material is given: $\int_\Omega\chi\,d\mathbf{x}=q_\alpha$. 
This problem could be understood
as the averaged version of the multiple state optimal design problem described 
above.

The case $l=1$ has been  studied in \cite{BM}. If $S$ consists of just one point, 
the problem reduces to  the (deterministic) multiple state optimal design problem 
which we started from. Therefore, the nonexistence of
a solution usually occurs, so the proper relaxation should be introduced.
Furthermore, in the case of a discrete probability space it is easily seen that 
an averaged multiple state optimal design problem can be
written as a deterministic multiple state optimal design problem, although 
the number of state equations significantly increases, even
to infinity in the case of an infinite probability space.

The content of the paper is as follows. Section 2 is devoted to description of 
the main issues regarding the proper relaxation
of the problem. The relaxed problem has a solution, which can be achieved as 
a (weak) limit of a minimizing sequence of the original problem, and vice versa: 
each minimizing sequence of the original problem has a subsequence converging to 
a solution of the relaxed problem \cite{MT}.
In Section 3 we further analyse the necessary conditions of optimality for the 
relaxed problem.

Section 4 deals with the special type of cost functional:
we consider the energy functional and maximise its average. 
An example is presented, where it is even possible to calculate the optimal design
explicitly by the technique presented in \cite{V2}. It shows that the averaged 
optimal control can differ significantly
from the optimal control obtained by averaged data, but this is due to 
the nonuniqueness of the averaged control.

\section{Proper relaxation}


If $\mathbf{A}\in\mathcal{M}(\alpha,\beta;\Omega)$ and $\mathcal{D} u:=-\operatorname{div}(\mathbf{A}\nabla u)$, then $\mathcal{D}$ is an isomorphism
between $ H_0^1(\Omega)$ and $H^{-1}(\Omega)$ by the Lax-Milgram lemma.
As an easy consequence of the following result
$\mathcal{D}$ is also an isomorphism between $L^p(S; H_0^1(\Omega))$ and 
$L^p(S;H^{-1}(\Omega))$, for any $p\in[1,+\infty]$.

\begin{lemma}\label{lplem}
If $X$ and $Y$ are two Banach spaces and $T$ is a continuous linear operator 
from $X$ to $Y$, then mapping $\mathcal{T}$ defined by
$$
(\mathcal{T}{\sf f})(s) := T({\sf f}(s))
$$
is a continuous linear operator from $L^p(S; X)$ to $L^p(S; Y)$.
\end{lemma}

The proof of Lemma \ref{lplem} follows from the construction of the 
B\^ochner integral: the case
$p=1$ is considered in \cite[Section V.5]{Y} or \cite[Section I.2]{DU}, 
while the general case can be treated in a similar vein.

For a sequence $(\mathbf{A}_n)$ in $\mathcal{M}(\alpha,\beta;\Omega)$, the corresponding sequence of isomorphisms
is introduced by $\mathcal{D}_n u:=-\operatorname{div}(\mathbf{A}_n\nabla u)$.
Following an analogous  approach for our original problem  we say
that a sequence $(\mathbf{A}_n)$ in $\mathcal{M}(\alpha,\beta;\Omega)$ converges  to $\mathbf{A}\in\mathcal{M}(\alpha',\beta';\Omega)$ if for
any $f\in L^p(S;H^{-1}(\Omega))$ and $u:=\mathcal{D}^{-1}f$, $u_n:=\mathcal{D}_n^{-1}f$  we have
\begin{gather*}
u_n \rightharpoonup u \quad \text{in } L^p(S; H_0^1(\Omega))\\
\mathbf{A}_n\nabla u_n \rightharpoonup \mathbf{A}  \nabla u \quad \text{in }
 L^p(S; L^2(\Omega;\mathbb{R}^d))\,.
\end{gather*}

The following Theorem establishes that this convergence actually coincides
 with classical $H$-convergence for stationary diffusion equation and with
 the notion of $H$-convergence  introduced in \cite{BM} 
(part 2 of Theorem \ref{th1}). The proof is straightforward.

\begin{theorem}\label{th1}
Let $\mathbf{A},\mathbf{A}_1,\mathbf{A}_2,\ldots$ be a sequence in $\mathcal{M}(\alpha,\beta;\Omega)$ and $p\in[1,+\infty\rangle$. 
The following statements are equivalent:
\begin{enumerate}
\item $\mathbf{A}_n$ $H$-converges to $\mathbf{A}$.

\item For any $f\in L^p(S;H^{-1}(\Omega))$ and $u:=\mathcal{D}^{-1}f$, $u_n:=\mathcal{D}_n^{-1}f$
 we have
\begin{gather*}
u_n(s,\cdot) \rightharpoonup u(s,\cdot) \quad \text{in }  H_0^1(\Omega)\\
\mathbf{A}_n\nabla u_n(s,\cdot) \rightharpoonup \mathbf{A}\nabla u(s,\cdot) \quad \text{in }
  L^2(\Omega;\mathbb{R}^d)\,,\quad {\rm for} \; {\rm a.e. }\; s\in S\,.
\end{gather*}

\item For any $f\in L^p(S;H^{-1}(\Omega))$ and $u:=\mathcal{D}^{-1}f$,
$u_n:=\mathcal{D}_n^{-1}f$  we have
\begin{gather*}
u_n \rightharpoonup u \quad \text{in } L^p(S; H_0^1(\Omega))\\
\mathbf{A}_n\nabla u_n \rightharpoonup \mathbf{A}  \nabla u \quad \text{in } 
L^p(S; L^2(\Omega;\mathbb{R}^d))\,.
\end{gather*}
\end{enumerate}
\end{theorem}


According to the theory of deterministic multiple state optimal design problems 
\cite{Ta2,A}, the  relaxation of the original
optimal design problem via homogenization theory consists in extending the
 original cost functional to
$$
J(\theta,\mathbf{A})=\int_S\int_\Omega \Bigl(\chi(\mathbf{x})\, g_\alpha(s,\mathbf{x},{\sf u}(s,\mathbf{x})) 
+ (1-\chi(\mathbf{x}))\, g_\beta(s,\mathbf{x},{\sf u}(s,\mathbf{x})\Bigr)\, d\mathbf{x}\,d\mu\,,
$$
where ${\sf u}$ is the unique solution of \eqref{ses}. The first step is to 
propose conditions which ensure the continuity of $J$ on $ L^\infty(\Omega)\times \mathcal{M}(\alpha,\beta;\Omega)$ 
in a reasonable topology.

To this end let us assume that ${\sf f}=(f_1,\ldots,f_m)\in L^p(S;H^{-1}(\Omega)^m)$ and
 $g_\alpha$ and $g_\beta$ to be Carath\'eodory functions (measurable in 
$(s,\mathbf{x})$ on product space $S\times\Omega$ 
and continuous in ${\sf u}$) satisfying the  growth condition:
\begin{equation}\label{growth1}
|g_\gamma(s,\mathbf{x},{\sf u})|\leq \varphi_\gamma(\mathbf{x},s) + \psi_\gamma(\mathbf{x},s){|{\sf u}|}^q 
\quad\text{for } \gamma=\alpha,\beta,
\end{equation}
with $\varphi_\gamma\in L^1(S\times\Omega)$, 
$\psi_\gamma\in L^{p'}(S; L^{q'}(\Omega))$, where $p'$ denotes the conjugate index 
to $p$: $\frac1{p}+\frac1{p'}=1$, and  $q\in [1,q^\ast\rangle$, with
$$
q^\ast= \begin{cases}
+\infty\,,& d\leq 2\\
{\frac{2d}{d-2}}\,, & d>2\,.
\end{cases}
$$

\begin{theorem}\label{thmcont}
Let the growth conditions \eqref{growth1} be satisfied for some $p\geq 1$. Then
for any ${\sf f}\in L^p(S;H^{-1}(\Omega) ^m)$
the functional $J$ is well defined and continuous on $ L^\infty(\Omega)\times \mathcal{M}(\alpha,\beta;\Omega)$
with respect to the weak-$\ast$ topology for $\theta$ and the $H$-topology 
for $\mathbf{A}$.
\end{theorem}

\begin{proof}
Let $f\in L^p(S;H^{-1}(\Omega)^m)$. If a sequence $(\mathbf{A}_n)$ in $\mathcal{M}(\alpha,\beta;\Omega)$ $H$-converges
to $\mathbf{A}\in\mathcal{M}(\alpha,\beta;\Omega)$, then by Theorem \ref{th1} for the sequence of
state functions ${\sf u}^i_n:=\mathcal{D}_n^{-1}f_i$ and ${\sf u}^i:=\mathcal{D}^{-1}f_i$,  $i=1,\ldots,m$, 
we have the weak convergence
\begin{equation}\label{slun}
u^i_n(s,\cdot)\rightharpoonup u^i(s,\cdot) \quad \text{in }
  H_0^1(\Omega)\,,\quad \text{for  a.e. } s\in S\,.
\end{equation}
By the  Sobolev imbedding theorem, the convergence is actually strong 
in $ L^q(\Omega)$, for any $q\in[1,q^\ast\rangle$ with $q^\ast$  given above.
Since $g_\alpha$ and $g_\beta$ are Carath\'eodory functions, this implies 
the convergence (up to a subsequence)
$$
g_\gamma (s,\mathbf{x},{\sf u}_n(s,\mathbf{x}))\to g_\gamma(s,\mathbf{x},{\sf u}(s,\mathbf{x})
$$
almost everywhere on $S\times\Omega$. By growth conditions \eqref{growth1} 
and the Lebesgue dominated convergence \cite[Section II.2]{DU} we
conclude that $g_\gamma (\cdot,\cdot,{\sf u}_n)$ converges to
 $g_\gamma (\cdot,\cdot,{\sf u})$ strongly in $ L^1(S\times\Omega)$, for 
$\gamma\in\{\alpha,\beta\}$.
Now if $\theta_n \smash{\rightharpoonup}^\ast\,
\theta$ in  $ L^\infty(\Omega)$ we conclude
 $J(\theta_n,\mathbf{A}_n)\to J(\theta,\mathbf{A})$.
Since the considered topology on $ L^\infty(\Omega)\times\mathcal{M}(\alpha,\beta;\Omega)$ is metrizable on
 bounded sets, this means that $J$ is continuous.
\end{proof}

The proper relaxation by homogenization theory consists in introducing 
the set of generalized designs
$$
\mathcal{A}:=\Big\{(\theta,\mathbf{A})\in{\rm L}^\infty(\Omega;[0,1]\times{\rm Sym_d}):
 \int_\Omega\theta\,d\mathbf{x}=q_\alpha,\,
\mathbf{A}(\mathbf{x})\in\mathcal{K}(\theta(\mathbf{x})) \text{ a.e. on }\Omega\Big\}\,,
$$
where $\mathcal{K}(\theta)$ stands for $G$-closure of the original set of conductivities 
$\mathbf{A}=\chi\alpha\mathbf{I}+(1-\chi)\beta\mathbf{I}$,
$\chi\in L^\infty(\Omega;\{0,1\})$, with given local fraction $\theta$ of the 
first phase.
The set $\mathcal{A}$ is compact with respect to the above considered topology,
and it is the closure of the set
of original designs. Therefore, the following result is straightforward.

\begin{theorem} \label{thm3}
The proper relaxation of the original problem reads
\begin{equation}\label{odprelax}
\begin{gathered}
J(\theta,\mathbf{A})=\int_S\int_\Omega \Bigl(\theta\, g_\alpha(\cdot,\cdot,{\sf u}) 
+ (1-\theta)\, g_\beta(\cdot,\cdot,{\sf u})\Bigr)\, d\mathbf{x}\,d\mu\to \min\\
(\theta,\mathbf{A})\in\mathcal{A}\,,\quad
{\sf u}=(u_1,\ldots, u_m)\text{ solves }\eqref{ses}\,.
\end{gathered}
\end{equation}
\end{theorem}

\begin{remark}\label{remar1} \rm
If we start from a conic sum of energies for each state equation and take 
its average over $S$:
$$
J(\theta,\mathbf{A})=\sum_{i=1}^m\lambda_i
 \int_S {}_{H^{-1}(\Omega)} \langle f(s),u_i(s)\rangle _{H_0^1(\Omega)} \,d\mu\,,
$$
it is sufficient to assume $f\in L^2(S;H^{-1}(\Omega))
=\bigl( L^2(S; H_0^1(\Omega))\bigr)'$
in order to obtain the continuity of $J$.
In particular, the relaxation problem has a form as written  in the previous Theorem.
\end{remark}

\section{Necessary conditions of optimality}

Let $({\theta^\ast},{\mathbf{A}^\ast})\in\mathcal{A}$ denote a solution of the relaxed problem \eqref{odprelax}
with corresponding state ${{\sf u}^\ast}$ and
let $\varepsilon\mapsto ({\theta^\varepsilon},{\mathbf{A}^\varepsilon})\in\mathcal{A}$ be a smooth path in $\mathcal{A}$ passing
through $({\theta^\ast},{\mathbf{A}^\ast})$ for $\varepsilon=0$.
By $u_i^\varepsilon$ we denote the corresponding state function, the unique solution of
\begin{equation}\label{seseps}
\begin{gathered}
-\operatorname{div}({\mathbf{A}^\varepsilon}\nabla u_i)=f_i(s,\cdot)\\
u_i \in L^p(S; H_0^1(\Omega))
\end{gathered}
\end{equation}
for $i=1\ldots m$.


After denoting  ${\delta\theta}=\frac{d}{d\varepsilon}{\theta^\varepsilon}\bigl|_{\varepsilon=0}$ and 
$\delta\mathbf{A}=\frac{d}{d\varepsilon}{\mathbf{A}^\varepsilon}\bigl|_{\varepsilon=0}$,
we would like  to calculate the variation of $J$ in terms of ${\delta\theta}$ and $\delta\mathbf{A}$, 
 more precisely, we look for
$\delta J:=\frac{d}{d\varepsilon}J({\theta^\varepsilon},{\mathbf{A}^\varepsilon})\bigl|_{\varepsilon=0}$:
\begin{equation}
\begin{aligned}
\delta J&= \lim_{\varepsilon\to0}\frac1\varepsilon
 \Big[\int_{S}\int_\Omega \theta^\varepsilon(\mathbf{x}) g_\alpha(s,\mathbf{x},{{\sf u}^\varepsilon}(s,\mathbf{x}))
 -{\theta^\ast}(\mathbf{x}) g_\alpha(s,\mathbf{x},{{\sf u}^\ast}(s,\mathbf{x}))\,d\mathbf{x}\,d\mu \\
&\quad + \int_{S}\int_\Omega(1-\theta^\varepsilon(\mathbf{x})) g_\beta(s,\mathbf{x},{{\sf u}^\varepsilon}(s,\mathbf{x}))
 -(1-{\theta^\ast}(\mathbf{x})) \\
&\quad\times g_\beta(s,\mathbf{x},{{\sf u}^\ast}(s,\mathbf{x}))\,d\mathbf{x}\,d\mu\Big]\\
&=\lim_{\varepsilon\to0}\Big[\int_{S}\int_\Omega
\frac{\theta^\varepsilon-{\theta^\ast}}\varepsilon (g_\alpha(\cdot,\cdot,{{\sf u}^\varepsilon})
 - g_\beta(\cdot,\cdot,{{\sf u}^\varepsilon})) \,d\mathbf{x}\,d\mu \\
&\quad + \int_{S}\int_\Omega {\theta^\ast}\sum_{i=1}^m
 \frac{u_i^\varepsilon-u_i^\ast}\varepsilon\int_0^1
 \frac{\partial g_\alpha}{\partial u_i}(\cdot,\cdot,{{\sf u}^\varepsilon_\tau})\,d\tau\,d\mathbf{x}\,d\mu\\
&\quad + \int_{S}\int_\Omega (1-{\theta^\ast})\sum_{i=1}^m
 \frac{u_i^\varepsilon-u_i^\ast}\varepsilon\int_0^1
 \frac{\partial g_\beta}{\partial u_i}(\cdot,\cdot,{{\sf u}^\varepsilon_\tau}))\,d\tau\,
\,d\mathbf{x}\,d\mu\Big]\,.
\end{aligned} \label{eta}
\end{equation}

Here, ${{\sf u}^\varepsilon_\tau}$ denotes ${{\sf u}^\ast}+\tau({{\sf u}^\varepsilon}-{{\sf u}^\ast})$, and it is assumed that 
$\frac{\partial g_\alpha}{\partial u_i}$ and $\frac{\partial g_\beta}{\partial u_i}$
are Carath\'eodory functions (measurable in $(s,\mathbf{x})$ and continuous in ${\sf u}$).
In order to pass to the limit above we additionally assume the following growth 
conditions
$$
\big|\frac{\partial g_\gamma}{\partial u_i}(s,\mathbf{x},{\sf u})\big|
\leq \widetilde\varphi_\gamma(s,\mathbf{x}) 
 + \widetilde\psi_\gamma(s,\mathbf{x}){|{\sf u}|}^r \quad{\rm for}\;
\gamma=\alpha,\beta;\; i=1,\ldots, m\,,
$$
where $\widetilde\varphi_\gamma\in L^{p'}(S; L^{q'}(\Omega))$,
$\widetilde\psi_\gamma\in L^{p_r}(S;L^{q_r}(\Omega))$, with   
$q$ being  the same as in the previous section,
$p_r=\frac{p}{p-r-1}$ and $q_r=\frac{q}{q-r-1}$
for some $r\geq 0$ such that $r\leq p-1$ and $r\leq q-1$.

By Theorem \ref{th1}, since L$^\infty$ convergence implies $H-$convergence, 
as in the proof of Theorem \ref{thmcont} we conclude that ${{\sf u}^\varepsilon}\to{{\sf u}^\ast}$
in $ L^q(\Omega)$ for any $q\in[1,q^\ast\rangle$, almost everywhere on 
$S$, as well as ${{\sf u}^\varepsilon_\tau}$, for any $\tau\in[0,1]$. Now, an application of the
 Lebesgue dominated convergence theorem implies that the
inner integrals over $\tau$ converge to 
$\int_0^1 \frac{\partial g_\gamma}{\partial u_i}(\cdot,\cdot,{{\sf u}^\ast})\,d\tau$
in $ L^{p'}(S;L^{q'}(\Omega))$, for $\gamma=\alpha,\beta$. If $\delta{\sf u}$ denotes 
$\frac{d}{d\varepsilon}{{\sf u}^\varepsilon}\bigl|_{\varepsilon=0}$, one concludes
\begin{align*}
{\delta J}&=\int_{S}\int_\Omega {\delta\theta}\bigl[g_\alpha(\cdot,\cdot,{{\sf u}^\ast})
 -g_\beta(\cdot,\cdot,{{\sf u}^\ast})\bigr] \\
&\quad + \sum_{i=1}^m \delta{\sf u}_i \Big({\theta^\ast}\frac{\partial g_\alpha}{\partial u_i}
 (\cdot,\cdot,{{\sf u}^\ast})+(1-{\theta^\ast})\frac{\partial g_\beta}{\partial u_i}
 (\cdot,\cdot,{{\sf u}^\ast})\Big)\,d\mathbf{x}\,d\mu \,. 
\end{align*}
Here, $\delta{\sf u}$ solves
\begin{gather*}
\operatorname{div}({\mathbf{A}^\ast}\nabla \delta u_i)
=\operatorname{div}(\delta\mathbf{A}\nabla u_i^\ast)\\
\delta u_i\in L^p(S; H_0^1(\Omega))
\end{gather*}
for $i=1,\ldots,m$.
We introduce adjoint states, as commonly, in order to eliminate those 
derivatives from the expression for ${\delta J}$.
Since $\theta\frac{\partial g_\gamma}{\partial u_i}$ belongs to 
$ L^{p'}(S;L^{q'}(\Omega))\hookrightarrow  L^{p'}(S;H^{-1}(\Omega))$, for
$\gamma=\alpha,\beta$,
the following boundary value problems have unique solutions
\begin{gather*}
-\operatorname{div}({\mathbf{A}^\ast}\nabla p_i)=\theta\frac{\partial g_\alpha}{\partial u_i}(\cdot,\cdot,{{\sf u}^\ast})
+(1-\theta)\frac{\partial g_\beta}{\partial u_i}(\cdot,\cdot,{{\sf u}^\ast})\\
p_i\in L^{p'}(S; H_0^1(\Omega))
\end{gather*}
for $i=1,\ldots, m$.
Now, one concludes that
\begin{equation}\label{gat}
{\delta J}=\int_{S}\int_\Omega {\delta\theta}\bigl[g_\alpha(\cdot,\cdot,{{\sf u}^\ast})
-g_\beta(\cdot,\cdot,{{\sf u}^\ast})\bigr]\,d\mathbf{x}\,d\mu
        -  \int_{S}\int_\Omega \sum_{i=1}^m\delta\mathbf{A}\nabla u_i^\ast\cdot
\nabla p_i^\ast\,d\mathbf{x}\,d\mu\,.
\end{equation}

\begin{remark}[Conic sum of energies] \rm
In the case of energy functional the calculation of variation ${\delta J}$ 
is straightforward, with the same assumption as in Remark \ref{remar1}:
$f\in L^2(S;H^{-1}(\Omega))$. The formula for the variation is the same,
 with ${p^\ast}={{\sf u}^\ast}$ in the case of minimization, and with ${p^\ast}=-{{\sf u}^\ast}$
for maximization problem.
\end{remark}

\begin{remark}[Discrete probability space] \rm
If $S$ is a finite set, the original random problem is easily seen as
 a (deterministic) multiple state optimal design problem.
The similar calculation holds for infinite, but discrete case: 
$S:=\{s_1,s_2,\ldots\}$, with probabilities $\mu(s_j)=\mu_j\geq 0$, 
$\sum_j\mu_j=1$.

Each $f_i\in L^p(S;H^{-1}(\Omega))$ determines the sequence of functionals
$f_i^j:=f(s_j,\cdot)\in H^{-1}(\Omega)$.
Denoting $u_i^j:=\mathcal{D}^{-1}(f_i^j)=(\mathcal{D}^{-1}f_i)(s_j,\cdot)$,
 ${\sf u}^j:=(u_1^j,\ldots u_m^j)$
and $g_\gamma^j(\mathbf{x},\sf{v})=g_\gamma(s_j,\mathbf{x},\sf{v})$, for 
$\gamma=\alpha,\beta$, we have
\begin{eqnarray*}
J(\theta,\mathbf{A}) 
&= \int_\Omega \theta\, \Big(\sum_j\mu_j g_\alpha^j(\cdot,{\sf u}^j)\Big) 
+ (1-\theta)\,\Big(\sum_j\mu_j
g_\beta^j(\cdot,{\sf u}^j)\Big)\, d\mathbf{x}\,.
\end{eqnarray*}
Here, the order of summation and integration can be interchanged by the
 Fubini theorem, due to assumptions \eqref{growth1}.

Denoting $\sf{U}=({\sf u}^1,{\sf u}^2,\ldots)$ and 
$h_\gamma(\mathbf{x},\sf{U})=\sum_j \mu_j g_\gamma^j(\mathbf{x},{\sf u}^j)$ (for $\gamma=\alpha, \beta$) 
we finally arrive at a multiple state optimal design problem 
(with infinite number of  state equations):
\begin{gather*}
J(\theta,\mathbf{A})=\int_\Omega \Bigl(\theta h_\alpha(\cdot,\sf{U}) 
+ (1-\theta) h_\beta(\cdot,\sf{U})\Bigr)\, d\mathbf{x}\to \min  \\
(\theta,\mathbf{A})\in\mathcal{A}\,.
\end{gather*}
\end{remark}

\section{Example}

Consider an energy maximization problem: take
$\Omega$ to be a ball $B(\mathbf{0},1)\subseteq\mathbb{R}^2$, $m=1$,
$q_\alpha:=0.8|\Omega|$, $S=\{1,2\}$ with $\mu_1=\mu_2=1/2$. 
The right-hand side is given by
$f_1:=f(1,\cdot)=\chi_{A} + \varepsilon\chi_B$ and
$f_2:=f(2,\cdot)=\chi_{A} - \varepsilon\chi_B$,
where
$A:=B(\mathbf{0},1/2)^c$ and $B:=B(\mathbf{0},1/5)$, as depicted in 
Figure \ref{primjer1}.

\begin{figure}[htb]
\begin{center}
 \includegraphics[scale=1]{fig1} % mv_ns151.mps
\end{center}
\caption{Right-hand sides are small perturbations of a constant heat 
source on annulus $B(\mathbf{0},1)\setminus B(\mathbf{0},1/2)$}
\label{primjer1}
\end{figure}

The average right-hand side $\overline f$ is simply $\chi_A$, and for a small 
$\varepsilon$ it would be interesting to compare solutions for the perturbed 
(right-hand sides are $f_1$ and $f_2$) and the unperturbed problem 
(right-hand side is $\overline f$).
By using the optimality criteria method \cite{A} we obtain numerical 
solutions presented in
Figure \ref{nums}.

\begin{figure}[htb]
\begin{center}
 \includegraphics[scale=0.2]{fig2a} % pr1_eps10-2.jpg
 \includegraphics[scale=0.2]{fig2b} % pr1_eps0.jpg
\end{center}
\caption{a. Numerical solution, $\varepsilon=0.01$;
b. Numerical solution, $\varepsilon=0$}\label{nums}
\label{fig2}
\end{figure}

For this example, due to spherical symmetry, it is possible to calculate the 
exact solution \cite{V2} for any $\varepsilon>0$ and examine what happens for
 $\varepsilon\to 0$. The method is based on the study of  necessary conditions of optimality \cite{MT}.
We shall just sketch the final result here.

Our problem (because of its spherical symmetry)
is equivalent to a simpler relaxation problem written only in
terms of local fraction $\theta$ \cite{V2}:
\begin{equation}\label{relm}
 I(\theta)
=\sum_{j=1}^2\mu_j\int_\Omega f_ju_j\,d\mathbf{x} \to \max 
\end{equation}
where $\theta\in L^\infty(\Omega;[0,1])$, $\int_\Omega\theta\,d\mathbf{x}=q_\alpha$,
and $u_1, u_2$ are determined uniquely by
\begin{gather*}
-\operatorname{div}(\lambda_-(\theta)\nabla u_j)=f_j\\
u_j\in {\rm H}^1_0(\Omega)
\end{gather*}
for $j=1,2$.
To be more precise,  for any solution $({\theta^\ast},{\mathbf{A}^\ast})$ of proper relaxation 
\eqref{odprelax}, ${\theta^\ast}$ is a solution of \eqref{relm}, and for any solution
$\widetilde{\theta}$ of \eqref{relm} we can construct a solution 
$\widetilde{\theta},\widetilde{\mathbf{A}}$ of \eqref{odprelax}
by taking simple laminate $\widetilde{\mathbf{A}}$ with local fraction $\widetilde{\theta}$ and 
layers orthogonal to the radial direction, at almost any point of 
domain \cite[Theorem 3.2]{V2}.

The necessary (and sufficient) conditions of optimality for \eqref{relm}
state  that there exists unique functions $\boldsymbol{\sigma}^*_i$
in $ L^2(\Omega;\mathbb{R}^2)$ satisfying $-\operatorname{div}\boldsymbol{\sigma}^*_i=f_i$ and
a Lagrange multiplier  $c\geq 0$ such that
\begin{equation}\label{nec2}
\begin{gathered}
\sum_{j=1}^2\mu_i|\boldsymbol{\sigma}^*_j|^2 > c \Rightarrow \theta^\ast=1\,,\\
\sum_{j=1}^2\mu_i|\boldsymbol{\sigma}^*_j|^2 < c \Rightarrow \theta^\ast=0\,.
\end{gathered}
\end{equation}
From the spherical symmetry one can show that $\boldsymbol{\sigma}^*_i$ are radial 
functions, and if $\Omega$ is a ball they can be uniquely determined by solving 
$-\operatorname{div}\boldsymbol{\sigma}^*_j=f_j$.

Here, we can calculate explicitly 
$\psi=\frac12(|\boldsymbol{\sigma}^*_1|^2+|\boldsymbol{\sigma}^*_2|^2)$
\[
\psi(r)=\begin{cases}
\varepsilon^2 r^2/4, &0\leq r\leq 1/5\\[4pt]
 \frac{ \varepsilon^2}{2500 r^2}, & 1/5< r \leq r 1/2\\[4pt]
 \frac{ r^2}4 + \big(\frac1{64}+\frac{\varepsilon^2}{2500}\big)\frac1{r^2}
-\frac18,  & 1/2 <r \leq 1
\end{cases}
\]
The graph of function $\psi$ is presented in Figure \ref{figpsi}. 
The optimality conditions \eqref{nec2} can  be used now to determine 
the unique optimal design,
as it is depicted in  Figure \ref{figpsi}.

\begin{figure}[htb]
\begin{center}
 \includegraphics[scale=1.2]{fig3} 
\end{center}
\caption{The graph of function $\psi=\frac12(|\boldsymbol{\sigma}^*_1|^2+|\boldsymbol{\sigma}^*_2|^2)$
and the geometric representation of optimality conditions.}
\label{figpsi}
\end{figure}

However, the limiting case $\varepsilon=0$ exhibits different behaviour: a solution
is not unique \cite{V2}: it is only important to set  $\alpha$ on 
$B(\mathbf{0},\frac12)^c$,
and to satisfy the constraint on the amount of the first phase.

\subsection*{Acknowledgments}
This work has been supported by Croatian Science Foundation under
the project 9780 WeConMApp.

The author wishes to thank the anonymous referee for the  comments 
and suggestions.


\begin{thebibliography}{99}

\bibitem{A} G.~Allaire; 
\emph{Shape optimization by the homogenization met\-hod}, Springer-Verlag, 2002.

\bibitem{AVmpe}  N.~Antoni\'{c}, M.~Vrdoljak; 
\emph{Sequential laminates in multiple state optimal design problems}, 
Math. Probl. Eng., {2006} (2006), Article ID 68695, 14 pages.

\bibitem{AVfer} N.~Antoni\'{c}, M.~Vrdoljak; 
\emph{Gradient methods for  multiple state optimal design problems},
Ann. Univ. Ferrara Sez. VII Sci. Mat., {53} (2007), pp. 177--187.

\bibitem{BM} G. Buttazzo, F. Maestre; 
\emph{Optimal shape for elliptic problems with random perturbations},
Discrete Contin. Dyn. Syst. 31 (2011), pp.~1115--1128.

\bibitem{DU} Joseph Diestel, John Jerry Uhl; 
\emph{Vector measures}, Amer. Math. Soc. 1977.

\bibitem{MT1} F.~Murat, L.~Tartar; 
\emph{H-convergence}, S\'eminaire d'Analyse Fonctionnelle et Num\'erique
 de l'Universit\'e d'Alger, mimeographed notes (1978), 
English translation in A. Cherkaev, R. V. Kohn (eds.); 
\emph{Topics in the mathematical modelling of composite materials}, 
Birkh\"auser, 1997.

\bibitem{MT} F.~Murat, L.~Tartar; 
\emph{Calcul des Variations et Homo\-g\'e\-n\'eisation}, in Les
 M\'ethodes de l'Homogenisation Th\'eorie et Applications
en Physique (Br\'eau-sans-Nappe, 1983), Collect. Dir. Etudes Rech. 
Elec. France {57},  Eyrolles, Paris, 1985, pp. 319--369.

\bibitem{Ssemel} S.~Spagnolo; 
\emph{Sul limite delle soluzioni di problemi di Cauchy relativi all'equazione 
del calore}, Ann. Scuola Norm. Sup. Pisa, {21}  (1967), pp. 657--699.

\bibitem{Ta1} L.~Tartar; 
\emph{Remarks on homogenization method in optimal design problems}, 
in Homogenization and applications to material sciences (Nice, 1995),
GAKUTO Internat. Ser. Math. Sci. Appl. {9}, Gakkokotosho, Tokyo, 1995, 
pp. 393--412.

\bibitem{Ta2} L.~Tartar; 
\emph{An introduction to the homogenization method in optimal design}, 
in Optimal shape design (Troia, 1998),
A.~Cellina and A.~Ornelas eds., {\it Lecture Notes in Math.}, 
\textbf{1740}, pp.~47--156, Springer-Verlag, 2000.

\bibitem{V} M.~Vrdoljak; 
\emph{On Hashin-Shtrikman bounds for mixtures of two isotropic materials},
 Nonlinear Anal. Real World Appl., {11} (2010),  pp. 4597--4606.

\bibitem{V2} M.~Vrdoljak; 
\emph{Classical optimal design in two-phase conductivity problems}, 
SIAM J. Control Optim., 54 (2016), pp. 2020--2035.

\bibitem{Y} K.~Yosida; 
\emph{Functional analysis}, Springer-Verlag, 1980.

\bibitem{ZKON} V. V.~\v{Z}ikov, S. M.~Kozlov, O. A.~Oleinik, H\`a Ti\^en Ngoan;
 \emph{Homogenization and G-convergence of differential operators},
Uspehi Mat. Nauk, {34} (1979), pp. 65--133.

\end{thebibliography}

\end{document}