Discretization approach

The discussion of discretization schemes is restricted to the standard situation where the domain is the unit square $\,\Omega=(0,1)\times (0,1)\,.$ The purpose of this section is to develop discretization techniques by which the distributed control problem (2.1)-(2.6) is transformed into a nonlinear programming problem (NLP-problem) of the form

$$\mbox{Minimize} \quad F^h(z) \quad \mbox{subject to} \quad G^h(z) = 0 \,, \quad H(z) \leq 0 \,.$$ (3.1)

The functions $\,F^h, G^h\,$ and $\,H\,$ are sufficiently smooth and are of appropriate dimension. The upper subscript $\,h\,$ denotes the dependence on the stepsize. The optimization variable $\,z\,$ will comprise both the discretized state and control variables.

The form (3.1) will be achieved by solving the elliptic equation (2.2) with the standard five-point-star discretization scheme. Choose a number $\,N \in I\!\!N_+\,$ and the stepsize $\,h:=1/(N+1)\,.$ Consider the mesh points

$$\,x_{ij}=(ih,jh)\,, \quad 0 \leq i,j \leq N+1, \,$$

and define the following sets of indices $\,(i,j)\,$ residing either in the domain $\,\Omega\,$ or on the boundary $\,\Gamma\,$, resp. on the subsets $\,\Gamma_1, \Gamma_2\,$ of the boundary:

$$\begin{array}{l} I(\Omega):=\{\,(i,j)\, \vert \; 1 \leq i,j ... ...mega\cup\Gamma_1):= I(\Omega) \cup I(\Gamma_1) \, . \end{array}$$ (3.2)

Obviously, we have $\;\char93 I(\Omega)= N^2\,, \;\char93 I(\Gamma)= 4*N\,$; define further $\;M_1:=\char93 I(\Gamma_1) \,$.

Now we shall present discretization schemes for the distributed control problem that are similar to those for boundary controls considered in Part 1 [23]. The optimization variable $\,z\,$ in (3.1) is taken as the vector

$$z:=(\,(y_{ij})_{\,(i,j) \in I(\Omega \cup \Gamma_1)}\,, \,(... ...j})_{\,(i,j) \in I(\Omega)}\,) \,\in I\!\! R^{\,2*N^2+M_1} \,.$$

The remaining state variables $\,y_{ij}, \, (i,j) \in I(\Gamma_2),$ are are determined by the Dirichlet condition (2.4) as

$$\hspace*{15mm} y_{ij} = y_2(x_{ij}) \quad \mbox{for} \quad (i,j) \in I(\Gamma_2) \,.$$ (3.3)

The derivative $\,\partial_{\nu} y(x_{ij}) \,$ in the direction of the outward normal is approximated by the expression $\, y^{\,\nu}_{ij}/h\,$ where

$$y^{\,\nu}_{ij} := \; \left \{ \begin{array}{llll} \ y_{i0} ... ... \mbox{for} & \; j=N+1, & i=1,...,N \end{array}\right \} \, .$$ (3.4)

Then the discrete form of the Neumann boundary condition (2.3) leads to the equality constraints

$$B^h(z):= y^{\nu}_{ij} - h\,b(x_{ij},y_{ij}) = 0 \qquad \mbox{for} \quad (i,j) \in I(\Gamma_1)\, . \quad$$ (3.5)

The application of the five-point-star to the elliptic equation $\, - \Delta y(x) + d(x,y(x),u(x)) \\ = 0 \,$ in (2.2) yields the following equality constraints for all $\,(i,j) \in I(\Omega)$:

$$\hspace*{-8mm} G^h_{ij}(z):= 4y_{ij} - y_{i+1,j} - y_{i-1,j}... ...i,j+1} - y_{i,j-1} + h^2\,d(x_{ij},y_{ij},u_{ij}) = 0 \,. \;$$ (3.6)

Note that the Dirichlet condition (3.3) is used in this equation to substitute the variables $\,y_{ij}\,$ for $\,(i,j) \in I(\Gamma_2)\,$. The control and state inequality constraints (2.5) and (2.6) yield the inequality constraints

$$C(x_{ij},y_{ij},u_{ij}) \leq 0 \,, \quad \forall \; (i,j) \in I(\Omega)\,,$$ (3.7)

$$S(x_{ij},y_{ij}) \leq 0 \,, \hspace*{12mm} \forall \; (i,j) \in I(\Omega\cup\Gamma_1)\,.$$ (3.8)

Observe that these inequality constraints do not depend on the meshsize $\,h\,$. The discretized form of the cost function (2.1) is

$$F^h(z):= h^2 \sum_{(i,j)\in I(\Omega)} f(x_{ij},y_{ij},u_{ij}) \,+\, h \sum_{(i,j) \in I(\Gamma_1)} g(x_{ij},y_{ij}) \,.$$ (3.9)

In summary, the relations (3.5)-(3.9) define a NLP-problem of the form (3.1). The corresponding Lagrangian function is

$$\hspace{-5mm} L(z,q,\lambda,\mu):= \hspace*{-3mm} h^2 \sum_{(i,j)\in I(\Omega)} f(x_{ij},y_{ij},u_{ij}) \,+\, h \sum_{(i,j) \in I(\Gamma_1)} g(x_{ij},y_{ij})$$ (3.10)

$$\hspace*{-25mm} +\,\sum_{(i,j)\in I(\Omega)} [\,q_{ij} G^h_{ij}(z) \, +\,\lambda_{ij} C(x_{ij},y_{ij},u_{ij}) + \,\mu_{ij} S(x_{ij},y_{ij})\,]$$

$$\hspace*{-25mm} +\, \sum_{(i,j)\in I(\Gamma_1)} [\,\mu_{ij} S(x_{ij},y_{ij})+q_{ij} B^h(z)\,],$$

where the Lagrange multipliers $\,q=(q_{ij})_{(i,j) \in I(\Omega\cup\Gamma_1)}\,$, $\,\lambda=(\lambda_{ij})_{(i,j) \in I(\Omega)}\,$ and $\,\mu=(\mu_{ij})_{(i,j) \in I(\Omega\cup\Gamma_1)}\,$ are associated with the equality constraints (3.5) and (3.6), resp. the inequality constraints (3.7) and (3.8). In addition, the multipliers $\,\lambda\,$ and $\,\mu\,$ satisfy the complementarity conditions corresponding to (2.14):

$$\begin{eqnarray*} \begin{array}{llrl} \lambda_{ij} \geq 0 & \quad\mbox{and} & \q... ... \quad \forall \; (i,j) \in I(\Omega\cup\Gamma_1) \,. \end{array}\end{eqnarray*}$$

In the next step we discuss the necessary conditions of optimality:

$$\,0 = L_z = (\,(L_{y_{ij}})_{(i,j)\in I(\Omega\cup\Gamma_1) }\,, (L_{u_{ij}})_{(i,j)\in I(\Omega)}\,) \,.$$

For state variables $\,y_{ij}\,$ with indices $\, (i,j) \in I(\Omega)\,$ we obtain the relations

$$\hspace*{-8mm} 0 = \hspace*{-4mm} L_{y_{ij}} = 4q_{ij} - q_{i+1,j} - q_{i-1,j} - q_{i,j+1} - q_{i,j-1} + h^2 q_{ij} \,d_y(x_{ij},y_{ij},u_{ij})\,+$$

$$\hspace*{3mm} +\,h^2 \,f_y(x_{ij},y_{ij},u_{ij}) \,+\, \lambda_{ij}\,C_y(x_{ij},y_{ij},u_{ij}) +\mu_{ij}\,S_y(x_{ij},y_{ij}) \,.$$ (3.11)

In these equations, the up to now undefined multipliers are set to

$$\hspace*{20mm} q_{ij} = 0 \qquad \forall \; (i,j) \in \Gamma_2 \,,$$ (3.12)

which is in accordance with the Dirichlet condition (2.12). We deduce from equations (3.11) that the Lagrange multipliers $\,q=(q_{ij})\,$ satisfy the five-point-star difference equations for the adjoint equation $\, -\Delta \bar{q} + \bar{q}\,d_y + f_y + \bar{\lambda}\,C_y + \,S_y\,\bar{\mu} = 0\,$ in (2.10) if we use the following approximations for the multiplier function $\,\bar{\lambda}\,$ and the regular Borel measure $\,\bar{\mu}\,$,

$$\hspace*{10mm} \bar{\lambda}(x_{ij}) \, \sim \, \lambda_{ij... ...\,, \quad \int_{sq(h^2)} \, d\bar{\mu} \; \sim \; \mu_{ij} \,,$$ (3.13)

where in the second relation $\,sq(h^2)\,$ denotes a square centered at $\,x_{ij}\,$ with area $h^2$. Recall the decomposition (2.15) of the measure $\,\bar{\mu} = \bar{\nu} \cdot dx\,+\,\bar{\nu}_s \cdot \bar{\mu}_s\,$. If the singular part of the measure vanishes, i.e. $\,\bar{\nu}_s \cdot \bar{\mu}_s=0\,$ holds, then (3.13) yields the following approximation for the density $\,\bar{\nu}$,

$$\hspace*{20mm} \bar{\nu}(x_{ij})\, \sim\ \mu_{ij}/h^2 \,.$$ (3.14)

In case that the measure $\, \bar{\mu} = \bar{\nu}_s \cdot \delta(x-x_{ij})\,$ is a delta distribution, we obtain from (3.13) the approximation

$$\hspace*{20mm} \bar{\nu}_s \, \sim\ \mu_{ij} \,.$$ (3.15)

For indices $\,(i,j) \in I(\Gamma_1)\,$ on the boundary $\,\Gamma_1\,$, e.g., for $\,j=0,\,i \in \{1,...,N\}\,$, we obtain

$$0= L_{y_{i0}} = -q_{i1} + q_{i0} -q_{i0}\,h\,b_y(x_{i0},y_{i0}) +h\,g_y(x_{i0},y_{i0}) + \mu_{ij}\,S_y(x_{ij},y_{ij}) \,.$$

These relations represent the discrete version of the Neumann boundary condition (2.11) if we approximate the regular Borel measure $\,\bar{\mu}\,$ on the boundary by

$$\hspace*{10mm} \int_{s(h)} \, d\bar{\mu} \; \sim \; \mu_{ij} \,,$$ (3.16)

where $\,s(h)\,$ denotes a line segment of length $\,h\,$ on $\,\Gamma_1\,$ centered at $\,x_{ij}\,.$ If the singular part in the decomposition (2.15) vanishes resp. if the measure $\, \bar{\mu} = \bar{\nu}_s \cdot \delta(x-x_{ij})\,$ is a delta distribution, we obtain the following approximations

$$\hspace*{10mm} \bar{\nu}(x_{ij})\, \sim\ \mu_{ij}/h \, \qquad \mbox{resp.} \quad \bar{\nu}_s \, \sim\ \mu_{ij} \,.$$ (3.17)

Finally, necessary conditions with respect to the control variables $\,u_{ij}\,$ for $\, (i,j) \in I(\Omega)\,$ are determined by

$$0= L_{u_{ij}} = h^2 f_u(x_{ij},y_{ij},u_{ij}) + q_{ij}\,h^2... ...,y_{ij},u_{ij}) + \lambda_{ij}\,C_u(x_{ij},y_{ij},u_{ij}) \,.$$

From this equation we recover the discrete version of the control law (2.13), if we use the same identification $\,\bar{\lambda}(x_{ij}) \,\sim\, \lambda_{ij}/h^2\,$ as in (3.13).