next up previous
Next: Optimization codes and modeling Up: Discretization and optimization techniques Previous: Discretization and optimization techniques

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 \,.
\end{displaymath} (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:
I(\Omega):=\{\,(i,j)\, \vert \; 1 \leq i,j ...
...mega\cup\Gamma_1):= I(\Omega) \cup I(\Gamma_1) \, .
\end{array}\end{displaymath} (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
y_{ij} = y_2(x_{ij}) \quad \mbox{for} \quad (i,j) \in I(\Gamma_2) \,.
\end{displaymath} (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 \{
\ y_{i0} ...
... \mbox{for} & \; j=N+1, & i=1,...,N
\end{array}\right \} \, .
\end{displaymath} (3.4)

Then the discrete form of the Neumann boundary condition (2.3) leads to the equality constraints
y^{\nu}_{ij} - h\,b(x_{ij},y_{ij}) = 0 \qquad \mbox{for} \quad
(i,j) \in I(\Gamma_1)\, . \quad
\end{displaymath} (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)$:
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 \,. \;
\end{displaymath} (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
    $\displaystyle C(x_{ij},y_{ij},u_{ij}) \leq 0 \,, \quad \forall \; (i,j) \in I(\Omega)\,,$ (3.7)
    $\displaystyle 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}) \,.
\end{displaymath} (3.9)

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

$\displaystyle \hspace{-5mm}
L(z,q,\lambda,\mu):=$   $\displaystyle \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)
    $\displaystyle \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})\,]$  
    $\displaystyle \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):

\lambda_{ij} \geq 0 & \quad\mbox{and} & \q...
... \quad \forall \; (i,j) \in I(\Omega\cup\Gamma_1) \,.

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
$\displaystyle \hspace*{-8mm}
0 =$   $\displaystyle \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})\,+$  
    $\displaystyle \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
q_{ij} = 0 \qquad \forall \; (i,j) \in \Gamma_2 \,,
\end{displaymath} (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}\,$,
\bar{\lambda}(x_{ij}) \, \sim \, \lambda_{ij...
...\,, \quad
\int_{sq(h^2)} \, d\bar{\mu} \; \sim \; \mu_{ij} \,,
\end{displaymath} (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}$,
\bar{\nu}(x_{ij})\, \sim\ \mu_{ij}/h^2 \,.
\end{displaymath} (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
\bar{\nu}_s \, \sim\ \mu_{ij}
\end{displaymath} (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
\int_{s(h)} \, d\bar{\mu} \; \sim \; \mu_{ij} \,,
\end{displaymath} (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
\bar{\nu}(x_{ij})\, \sim\ \mu_{ij}/h \,
\qquad \mbox{resp.} \quad \bar{\nu}_s \, \sim\ \mu_{ij} \,.
\end{displaymath} (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}) +
...,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).
next up previous
Next: Optimization codes and modeling Up: Discretization and optimization techniques Previous: Discretization and optimization techniques
Hans D. Mittelmann