# Numerical Shape Optimization via Dynamic Programming

Jan Pustelnik

► To cite this version:

Jan Pustelnik. Numerical Shape Optimization via Dynamic Programming. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.417-426, 10.1007/978-3-642-36062-6_42 . hal-01347564

HAL Id: hal-01347564

https://inria.hal.science/hal-01347564

Submitted on 21 Jul 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical shape optimization via dynamic programming

Jan Pustelnik

University of Lodz, Fac. of Math. & Computer Science,
Banacha 22, 90-128 Lodz, Poland

**Abstract.** In this paper we describe a novel framework for finding numerical solutions to a wide range of shape optimization problems. It is based on classical dynamic programming approach augmented with discretization of the space of trajectories and controls. This allows for straightforward algorithmic implementation. This method has been used to solve a well known problem called the "dividing tube problem", a state problem related to fluid mechanics, that requires simultaneous topology and shape optimization in case of elastic contact problems and involves solving the Navier-Stokes equations for viscous incompressible fluids.

**Keywords:** dynamic programming, numerical approximation, contact problem, shape optimization, sufficiency optimality condition, structural optimization, topological derivative, stationary Navier-Stokes equations.

## 1 Introduction

In the paper, as a model problem, we consider state problems related to fluid mechanics, namely the Navier-Stokes equations for viscous incompressible fluids. The main problem is to search for optimal shape of a given objective. For an incompressible fluid, conservation laws for momentum and mass are assumed to be in force. The displacement field of the body is governed by the Reynolds-averaged Navier-Stokes (RANS) equations with an algebraic mixing length turbulence. The volume of the body is assumed to be bounded.

The results pertaining the existence, regularity and finite-dimensional approximation of solutions to mentioned problems are given in [4], [5]. The primal-dual algorithms for numerical solving of contact problems were developed in [6], [8]. In the course of solution the necessary optimality condition for simultaneous shape and topology optimization is formulated, while the shape and topological derivatives are employed, what stays close to classical optimization problems and gives sufficient optimality conditions. It is a different approach than the one applied e.g. in [2] where the notion of topological derivative and results concerning its application in optimization of elastic structures are reported.

We describe a new numerical algorithm for that optimization problem.
## 2 General shape optimization problem

We consider the following shape optimization problem, which is being analyzed and subsequently solved in [3]:

$$ \text{minimize} \quad J(\Omega) = \int_{\Omega} L(x, u(x), \nabla u(x))dx \quad (1) $$

subject to

$$ \Omega \in \Theta, \quad Au(x) = f(x, u(x)), \quad Bu(x) = \phi(x) \quad \text{on} \quad \partial\Omega \quad (2) $$

where $\Theta$ is a certain family of bounded with $C^{0,1}$ boundary domains of $D \subset R^n$ which will be defined precisely in subsection 2.1 and $A$ is a differential operator e.g. defining Navier-Stokes equations and $B$ an operator acting on the boundary. We assume that $L: R^n \times R \times R^n \to R$ is Lipschitz continuous with respect to all variables, $f: R^n \times R^m \to R^m$ is continuous and Lipschitz continuous with respect to last variable, $\phi(\cdot)$ is continuous on $\partial\Omega$.

### 2.1 Reduction of shape optimization problem to classical control problem

We will summarize in this subchapter the results of research published in [3] in order to introduce relevant objects on which the presented numerical method operates.

Let $U$ be a given nonempty, compact set in $C^{0,1}$ of surfaces defined on $E \subset R^{n-1}$. We assume that each supremum of each subfamily of $U$ also belongs to $U$ as well as finite concatenation of element of $U$ belongs to $U$. Let $U \ni v \to \Omega(v)$ be a given family of simply connected domains in $D \subset R^n$ with $C^{0,1}$ boundary such that some fixed part of their boundary is changing and is a surface $v$ from $U$. We assume that $\Omega(v)$ depends in a smooth way on $v$ and that there exists a $v_{\max} \in U$ such that $\Omega(v) \subset \Omega(v_{\max})$, for all $v \in U$ and there exists a $v_{\min} \in U$ such that $\Omega(v_{\min}) \subset \Omega(v)$, for all $v \in U$. Let us denote the part of the boundary $\partial\Omega(v_{\min})$ corresponding to the surface $v_{\min}$ as $\Gamma_0$ while that corresponding the surface $v$ as $\Gamma(v)$. We have the following BVP:

Find $z_{\max} \in C^{1,k}(\Omega(v_{\max}))$ such that $\Delta z(x) = 0$ in $\Omega(v_{\max}) \setminus \bar{\Omega}(v_{\min})$, $z(x) = 0$ on $\Gamma_0$, $z(x) = 1$ on $\Gamma(v_{\max})$. Next for $v \in U$, $v \neq v_{\min}$, find $z \in C^{1,k}(\Omega(v) \setminus \bar{\Omega}(v_{\min}))$ such that: $\Delta z(x) = 0$ in $\Omega(v_{\max}) \setminus \bar{\Omega}(v_{\min})$, $z(x) = 0$ on $\Gamma_0$, $z(x) = z_{\max}(x)$ on $\Gamma(v)$. Solutions to the BVP above belong to $C^{1,k}(\Omega(v) \setminus \bar{\Omega}(v_{\min}))$ and in fact they are restrictions of $z_{\max}$ to $\Omega(v) \setminus \bar{\Omega}(v_{\min})$. Because $z(x)$ depends (in a continuous way) on $v$, we will use the notation $z(x, v)$. We define the family $\Theta$ of sets over which the problem (1)-(2) is considered as: $\Theta = \{\Omega(v) : v \in U\}$. The sets from $\Theta$ are called admissible sets for problem (1)-(2). For a given $\Omega(v) \setminus \bar{\Omega}(v_{\min})$, we introduce the field and the deformation: $V(x, v) = ||\nabla z(x, v)||^{-2}\nabla z(x, v)$, $T(w, v) = x(s, w, v)$, $s \in [0, 1]$, where $x(\cdot, w, v)$ is a solution to $\frac{d}{ds}x(s, w, v) = V(x(s, w, v), v)$, $s \in [0, 1]$, $x(0, w, v) = w$,
$w \in \Gamma_0$. Notice, that for a given fixed $w \in \Gamma_0$, the point $x(1, w, v)$ belongs to $\Gamma(v)$.

For a given control $v \in U$ we can write:

$$ \frac{d}{ds}x(s, w, v) = V(x(s, w, v), v), \quad s \in [0, 1], \quad x(0, w, v) = w. \qquad (3) $$

Then the boundary $\Gamma(v)$ is the image of $\Gamma_0$ by the map $x(1, \cdot, v)$. Thus, for a given $v \neq v_{\min}$, we have an alternative definition of the $\Omega(v) \setminus \bar{\Omega}(v_{\min})$: $\Omega(v) \setminus \bar{\Omega}(v_{\min}) = \{x: x = x(s, w, v), 0 < s < 1, w \in \Gamma_0\}$. This means that we can construct and study some objects over the set $\Omega(v)$ with the help of the family $F(v): F(v) = \{x(s, w, v): 0 < s < 1, w \in \Gamma_0\}$. The original functional $J(\Omega)$ in terms of the family $F(v)$ can be rewritten as $J(F(v)) = I(v)$,

$$ \mathbf{I}(v) = \int_{\Omega(v_{\min})} L(y, u(y), \nabla u(y))dy + \int_{\Omega(v) \setminus \bar{\Omega}(v_{\min})} L(x, u(x), \nabla u(x))dx \\ = \int_{\Omega(v_{\min})} L(y, u(y), \nabla u(y))dy + \int_{0}^{1} \int_{\Gamma_0} \hat{L}(x(s, w, v))dwds, $$

where $\hat{L}(x(s, w, v)) = L(x(s, w, v), u(x(s, w, v)), \nabla u(x(s, w, v))) \left| \frac{\partial}{\partial s}x \frac{\partial}{\partial w}x \right|$.

Therefore we are able to reduce the original shape optimal control problem to classical optimal control problem (P): minimize $\mathbf{I}(v)$ subject to $\frac{d}{ds}x(s, w, v) = V(x(s, w, v), v), s \in [0, 1], x(0, w, v) = w, w \in \Gamma_0, v \in U$,

$$ \Omega(v) \in \Theta, \quad \mathbb{A}u(x) = f(x, u(x)), \quad x \in \Omega(v). \qquad (4) $$

In order to formulate any sufficient optimality conditions for this problem we apply classical dynamic programming scheme.

## 2.2 Dynamic programming approach as a method to solution of (P)

Let us take any $x \in \Omega(v_{\max}) \setminus \bar{\Omega}(v_{\min})$ and denote by $U_x$ a subfamily of $U$ such that $x \in v$ for each $v \in U_x$. Next denote by $v_x = \max U_x$, where the maximum over $U_x$ means that $\Omega(v) \subset \Omega(v_x)$ for all $v \in U_x$. By our assumption on $U$, $v_x$ exists and $v_x \in U$. Put $\bar{U}_x = \{v \in U: \Omega(v) \subset \Omega(v_x)\}$. By (3) for each $v \in U_x$ there is a trajectory $x(\cdot, w, v)$ such that $x = x(1, w_v, v)$, for some $w_v \in \Gamma_0$. The problem (P) falls into the category of Lagrange control problems treated in many books (e.g. [1]). Following Chapter IV of this book we define a value function for (P), for $x \in \Omega(v_{\max})$:

$$ S(x) = \inf \left\{ \int_{\Omega(v_{\min})} L(y, u(y), \nabla u(y))dy + \int_{0}^{1} \int_{\Gamma_0} \hat{L}(x(s, w, v))dwds \right\}, \quad (5) $$

where infimum in (5) is taken over all pairs $(x(\cdot, w, v), v)$ satisfying $\frac{d}{ds}x(s, w, v) = V(x(s, w, v), v), s \in [0, 1], v \in \bar{U}_x, w \in \Gamma_0$ and for $v \in U_x, x(1, w_v, v) = x$, for some $w_v \in \Gamma_0$. Each pair $(x(\cdot, w, v), v)$ satisfying these equations will be called
admissible for the point $x \in \Omega(v_{\max})\setminus\bar{\Omega}(v_{\min})$. However, in practice, we cannot expect that $S(\cdot)$ is of $C^1$ in $\Omega(v_{\max})\setminus\bar{\Omega}(v_{\min})$, this is why we are interested in numerical approximation of $S(\cdot)$. Therefore, we shall look for $\varepsilon$-value function: $S_{\varepsilon}(\cdot)$. For given $\varepsilon > 0$ we call $S_{\varepsilon}: \Omega(v_{\max})\setminus\bar{\Omega}(v_{\min}) \to R$, $\varepsilon$-value function if

$$S(x) \leq S_{\varepsilon}(x) \leq S(x) + \varepsilon, x \in \Omega(v_{\max})\setminus\bar{\Omega}(v_{\min}). \quad (6)$$

It is clear that there exists infinitely many $\varepsilon$-value functions $S_{\varepsilon}(\cdot)$.

### 3 Numerical approximation of the value function

This section is an adaptation of the method developed by Pustelnik in his Ph.D. thesis [7] for numerical approximation of value function for Bolza problem from optimal control theory.

Let us define the following set $T = \{x: x \in \Omega(v_{\max})\setminus\bar{\Omega}(v_{\min})\}$. Since $\Omega(v_{\max})\setminus\bar{\Omega}(v_{\min})$ is bounded, the set $\bar{T}$ is compact. Let $T \ni x \to g(x)$ be an arbitrary function of class $C^1$ in $\bar{T}$ such that $g(x) = c, x \in \Gamma_0$, where $c$ is some constant which will be determined later. For a given function $g$, we define $(x, v) \to G_g(x, v)$ as

$$G_g(x, v) = g_x(x)V(x, v) - \int_{\Gamma_0} \hat{L}(x(1, w, v))dw, \quad (7)$$

$v \in \bar{U}_x$, where $x(\cdot, w, v)$, $u$ are defined as previously. Next, we define the function $x \to F_g(x)$ as

$$F_g(x) = \max \{G_g(x, v) : v \in \bar{U}_x\}. \quad (8)$$

Note that by the assumptions on $L$ and $V$, the function $F_g$ is continuous in $T$. By the continuity of $F_g$ and compactness of $\bar{T}$, there exist $k_d$ and $k_g$ such that $k_d \leq F_g(x) \leq k_g$ for $x \in \Omega(v_{\max})\setminus\bar{\Omega}(v_{\min})$.

#### 3.1 Definition of covering of $T$

Let $\eta > 0$ be fixed and $\{q_j^\eta\}_{j \in \mathbb{Z}}$ be a sequence of real numbers such that $q_j^\eta = j\eta$, $j \in \mathbb{Z}$ ($\mathbb{Z}$ - set of integers). Denote

$$J = \{j \in \mathbb{Z}: \text{ there is } x \in T, j\eta < F_g(x) \leq (j+1)\eta\},$$

Next, let us divide the set $T$ into the sets $P_j^{\eta,g}$, $j \in J$, as follows

$$P_j^{\eta,g} := \{x \in T : q_j^\eta < F_g(x) \leq q_{j+1}^\eta\}, \quad j \in J.$$

#### 3.2 Discretization of $F_g$

Define in $T$ a function

$$h^{\eta,g}(x) = -q_{j+1}^{\eta}, \quad x \in P_j^{\eta,g}, \quad j \in J. \quad (9)$$
Then, by the construction of the covering of T, we have

$$0 \le F_g(x) + h^{\eta,g}(x) \le \eta, \quad x \in T. \tag{10}$$

Let $(x(\cdot, w, v), v)$ be any admissible pair with the trajectory defined in $[0, 1]$, starting at the point $x(0, w, v)$, $w \in \Gamma_0$ fixed. We show that there exists an increasing sequence of $m$ points $\{\tau_i\}_{i=1,...,m}$, $\tau_1 = 0$, $\tau_m = 1$, such that for $\tau \in [\tau_i, \tau_{i+1}]$

$$|F_g(x(\tau_i, w, v)) - F_g(x(\tau, w, v))| \le \frac{\eta}{2}, \quad i = 2, \dots, m-2, \tag{11}$$

$$|F_g(x(\tau_2, w, v)) - F_g(x(\tau, w, v))| \le \frac{\eta}{2}, \quad \tau \in (\tau_1, \tau_2],$$

$$|F_g(x(\tau_{m-1}, w, v)) - F_g(x(\tau, w, v))| \le \frac{\eta}{2}, \quad \tau \in [\tau_{m-1}, \tau_m).$$

Indeed, it is a direct consequence of two facts: Lipschitz continuity of $x(\cdot, w, v)$ with a common Lipschitz constant and continuity of $F_g$. From (11) we infer that for each $i \in \{1, \dots, m-1\}$ if $x(\tau_i, w, v) \in P_j^{\eta,g}$ for a certain $j \in J$, then we have for $\tau \in [\tau_i, \tau_{i+1})$

$$x(\tau, w, v) \in P_{j-1}^{\eta,g} \cup P_j^{\eta,g} \cup P_{j+1}^{\eta,g}.$$

Define

$$h^{\eta,g}(x(\tau_1, w, v)) = h^{\eta,g}(x(\tau, w, v)) \quad \text{for some } \tau \text{ near } \tau_1,$$

$$h^{\eta,g}(x(\tau_m, w, v)) = h^{\eta,g}(x(\tau, w, v)) \quad \text{for some } \tau \text{ near } \tau_m.$$

Thus for $\tau \in [\tau_i, \tau_{i+1}]$

$$h^{\eta,g}(x(\tau_i, w, v)) - \eta \le h^{\eta,g}(x(\tau, w, v)) \le h^{\eta,g}(x(\tau_i, w, v)) + \eta, \tag{12}$$

and so, for $i \in \{2, \dots, m-1\}$

$$h^{\eta,g}(x(\tau_i, w, v)) - h^{\eta,g}(x(\tau_{i-1}, w, v)) = \eta_i^i(x(\cdot, w, v)), \tag{13}$$

where $\eta_i^i(x(\cdot, w, v))$ is equal to $-\eta$ or 0 or $\eta$. Integrating (12) we get

$$\left| \int_0^1 h^{\eta,g}(x(\tau, w, v)) d\tau - \sum_{i \in \{1, \dots, m-1\}} \left[ h^{\eta,g}(x(\tau_i, w, v))(\tau_{i+1} - \tau_i) \right] \right| \le \eta.$$

By using the formula above and the following simple arithmetic transformations

$$\sum_{i \in \{2, \dots, m-1\}} [h^{\eta,g}(x(\tau_i, w, v)) - h^{\eta,g}(x(\tau_{i-1}, w, v))](\tau_m - \tau_i) = \sum_{i \in \{1, \dots, m-1\}} [h^{\eta,g}(x(\tau_i, w, v))(\tau_{i+1} - \tau_i)] - h^{\eta,g}(x(\tau_1, w, v))(\tau_1 - \tau_m),$$
we obtain

$$
\begin{align*}
& \sum_{i \in \{2, \dots, m-1\}} [h^{\eta,g}(x(\tau_i, w, v)) - h^{\eta,g}(x(\tau_{i-1}, w, v))](\tau_m - \tau_i) \\
& + h^{\eta,g}(x(\tau_1, w, v))(\tau_m - \tau_1) - \eta(\tau_m - \tau_1) \\
& \le \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w, v)) d\tau \\
& \le \sum_{i \in \{2, \dots, m-1\}} [h^{\eta,g}(x(\tau_i, w, v)) - h^{\eta,g}(x(\tau_{i-1}, w, v))](\tau_m - \tau_i) \\
& + h^{\eta,g}(x(\tau_1, w, v))(\tau_m - \tau_1) + \eta(\tau_m - \tau_1).
\end{align*}
$$

and, taking into account (13), we infer that

$$
\begin{align}
& \sum_{i \in \{2, \dots, m-1\}} \eta_{x(\cdot, w, v)}^i (\tau_m - \tau_i) + h^{\eta,g}(x(\tau_1, w, v))(\tau_m - \tau_1) - \eta(\tau_m - \tau_1) \nonumber \\
& \le \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w, v)) d\tau \tag{14} \\
& \le \sum_{i \in \{2, \dots, m-1\}} \eta_{x(\cdot, w, v)}^i (\tau_m - \tau_i) + h^{\eta,g}(x(\tau_1, w, v))(\tau_m - \tau_1) + \eta(\tau_m - \tau_1). \nonumber
\end{align}
$$

We would like to stress that (14) is very useful from numerical point of view: we can estimate the integral $h^{\eta,g}(\cdot, \cdot)$ along any trajectory $x(\cdot, w, v)$ as a sum of finite number of values, where each value consists of a number from the set $\{-\eta, 0, \eta\}$ multiplied by $\tau_m - \tau_i$. Moreover, for two different trajectories: $x(\cdot, w^1, v^1)$, $x(\cdot, w^2, v^2)$, the expressions

$$
\sum_{i \in \{2, \dots, m-1\}} \eta_{x(\cdot, w^1, v^1)}^i (\tau_m - \tau_i) + h^{\eta,g}(x(\tau_1, w^1, v^1))(\tau_m - \tau_1)
$$

and

$$
\sum_{i \in \{2, \dots, m-1\}} \eta_{x(\cdot, w^2, v^2)}^i (\tau_m - \tau_i) + h^{\eta,g}(x(\tau_1, w^2, v^2))(\tau_m - \tau_1)
$$

are identical if

$$
h^{\eta,g}(x(\tau_1, w^1, v^1)) = h^{\eta,g}(x(\tau_1, w^2, v^2)) \quad (15)
$$

and

$$
\eta_{x(\cdot, w^1, v^1)}^{i} = \eta_{x(\cdot, w^2, v^2)}^{i} \text{ for all } i \in \{2, \dots, m-1\}. \quad (16)
$$

The last one means that in the set B of all trajectories $x(\cdot, w, v)$, $w \in \Gamma_0$, $v \in U$, we can introduce an equivalence relation $r$: we say that two trajectories $x(\cdot, w^1, v^1)$ and $x(\cdot, w^2, v^2)$, $w^1, w^2 \in \Gamma_0$, $v^1, v^2 \in U$ are equivalent if they satisfy (15) and (16). We denote the set of all disjoint equivalence classes by $B_r$. The cardinality of $B_r$, denoted by $||B_r||$, is finite and bounded from above by $3^{m+1}$.

Define

$$
X = \{ x = (x_1, \dots, x_{m-1}) : x_1 = 0, x_i = \eta_{x_j}^i, \\
\qquad i = 2, \dots, m-1, x^j \in B_r, j = 1, \dots, ||B_r|| \}.
$$
It is easy to see that the cardinality of X is finite.

The considerations above allow us to estimate the approximation of the value function.

**Theorem 1.** We have the following estimation

$$
\begin{aligned}
& \min_{x \in B_r, w_0 \in \Gamma_0} \left( - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, v)) d\tau - g(x(\tau_m, w_0, v)) \right) \\
& \leq \max_{x \in B_r} \left\{ \int_{\tau_1}^{\tau_m} \left( - \int_{\Gamma_0} \hat{L}(x(s, w, v)) dw \right) ds - g(x(\tau_1, w_0, v)) \right\} \\
& \leq \max_{x \in B_r, w_0 \in \Gamma_0} \left( - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, v)) d\tau - g(x(\tau_m, w_0, v)) \right) + \eta(\tau_m - \tau_1),
\end{aligned}
$$

where u is a solution to (4) for $\Omega(v)$.

*Proof.* By inequality (10) $0 \leq F_g(x) + h^{\eta,g}(x) \leq \eta$ we have $-h^{\eta,g}(x) \leq F_g(x) \leq -h^{\eta,g}(x) + \eta$. Integrating the last inequality along any $x(\cdot, w_0, \tilde{v})$ in the interval $[\tau_1, \tau_m]$ we get

$$
\begin{aligned}
& - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, \tilde{v})) d\tau \\
\leq & \max_{v \in U} \int_{\tau_1}^{\tau_m} \left( g_x(x(\tau, w_0, \tilde{v})) V(x(\tau, w_0, \tilde{v}), v) - \int_{\Gamma_0} \hat{L}(x(\tau, w, v)) dw \right) d\tau \\
\leq & - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, \tilde{v})) d\tau + \eta(\tau_m - \tau_1).
\end{aligned}
$$

Hence, we get two inequalities

$$
\begin{aligned}
& \min_{x \in B_r, w_0 \in \Gamma_0} \left( - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, \tilde{v})) d\tau - g(x(\tau_m, w_0, \tilde{v})) \right) \\
& \leq \min_{x \in B_r, w_0 \in \Gamma_0} \max_{v \in U} \int_{\tau_1}^{\tau_m} \left( -g(x(\tau_m, w_0, \tilde{v})) \right. \\
& \qquad \left. + g_x(x(\tau, w_0, \tilde{v})) V(x(\tau, w_0, \tilde{v}), v) - \int_{\Gamma_0} \hat{L}(x(\tau, w, v)) dw \right) d\tau
\end{aligned}
$$

and

$$
\begin{aligned}
& \max_{x \in B_r, w_0 \in \Gamma_0} \max_{v \in U} \int_{\tau_1}^{\tau_m} \left( -g(x(\tau_m, w_0, \tilde{v})) \right. \\
& \qquad \left. + g_x(x(\tau, w_0, \tilde{v})) V(x(\tau, w_0, \tilde{v}), v) - \int_{\Gamma_0} \hat{L}(x(\tau, w, v)) dw \right) d\tau \\
\leq & \max_{x \in B_r, w_0 \in \Gamma_0} \left( - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, \tilde{v})) d\tau - g(x(\tau_m, w_0, \tilde{v})) \right) + \eta(\tau_m - \tau_1).
\end{aligned}
$$
As a consequence of the above we get

$$
\begin{aligned}
& \min_{x \in B_r, w_0 \in \Gamma_0} \left( - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, \tilde{v})) d\tau - g(x(\tau_m, w_0, \tilde{v})) \right) \\
& \le \max_{x \in B_r} \left\{ \int_{\tau_1}^{\tau_m} \left( - \int_{\Gamma_0} \hat{L}(x(\tau, w, v)) dw \right) d\tau - g(x(\tau_1, w_0, v)) \right\} \\
& \le \max_{x \in B_r, w_0 \in \Gamma_0} \left( - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, v)) d\tau - g(x(\tau_m, w_0, \tilde{v})) \right) + \eta(\tau_m - \tau_1)
\end{aligned}
$$

and thus the assertion of the theorem follows.

Let us now define four following symbols: $\mathcal{F}_{(\eta,1)}(x) := - \sum_{i=2,...,m-1} \eta_i^i(\tau_m - \tau_1)$, $\mathcal{F}_1(x) := - \sum_{i \in \{1,...,m-1\}} x^i(\tau_m - \tau_1)$, $\mathcal{F}_{(\eta,i)}(x) := - \sum_{i=2,...,m-1} \eta_i^i(\tau_m - \tau_i)$, $\mathcal{F}_i(x) := - \sum_{i \in \{1,...,m-1\}} x^i(\tau_m - \tau_i)$. Now, we use the definition of equivalence class to reformulate the theorem above in a way that is more useful in practice. To this effect let us note that, by definition of equivalence relation $r$, we have

$$ \min_{x \in B_r} \{ \mathcal{F}_{(\eta,1)}(x) \} = \min_{x \in X} \{ \mathcal{F}_1(x) \}, \max_{x \in B_r} \{ \mathcal{F}_{(\eta,1)}(x) \} = \max_{x \in X} \{ \mathcal{F}_1(x) \}. $$

Let us now also define the following auxiliary symbol $\mathcal{H}(x, w_0) := -h^{\eta,g}(x(\tau_1, w_0, v))(\tau_m - \tau_1) - g(x(\tau_m, w_0, v))$. Taking into account (14) we get

$$
\begin{aligned}
& \min_{x \in X} \{ \mathcal{F}_i(x) \} + \min_{x \in B_r, w_0 \in \Gamma_0} \{ \mathcal{H}(x, w_0) \} - \eta(\tau_m - \tau_1) \\
& \le \min_{x \in B_r} \left\{ - \int_{\tau_1}^{\tau_m} h^{\eta,g}(x(\tau, w_0, v)) d\tau - g(x(\tau_m, w_0, v)) \right\} \\
& \le \min_{x \in X} \{ \mathcal{F}_i(x) \} + \max_{x \in B_r, w_0 \in \Gamma_0} \{ \mathcal{H}(x, w_0) \} + \eta(\tau_m - \tau_1)
\end{aligned}
$$

and a similar formula for supremum. Applying that to the result of the theorem above, we obtain the following estimation

$$
\begin{aligned}
& \min_{x \in X} \{ \mathcal{F}_i(x) \} + \min_{x \in B_r, w_0 \in \Gamma_0} \{ \mathcal{H}(x, w_0) \} - 2\eta(\tau_m - \tau_1) \\
& \le \max_{x \in B_r} \left\{ \int_{\tau_1}^{\tau_m} \left( - \int_{\Gamma_0} \hat{L}(x(\tau, w, v)) dw \right) d\tau - g(x(\tau_1, w_0, v)) \right\} \quad (17) \\
& \le \max_{x \in X} \{ \mathcal{F}_i(x) \} + \max_{x \in B_r, w_0 \in \Gamma_0} \{ \mathcal{H}(x, w_0) \} + \eta(\tau_m - \tau_1).
\end{aligned}
$$

Thus, we come to the main theorem of this section, which allows us to reduce an infinite dimensional problem to the finite dimensional one.

**Theorem 2.** Let $\eta > 0$. Assume that there is $\theta > 0$ and $\bar{v}$ such that

$$
\begin{aligned}
& \max_{x \in X} \{ \mathcal{F}_i(x) \} \\
& \quad + \max_{x \in B_r, w_0 \in \Gamma_0} \{ \mathcal{H}(x, w_0) \} \\
& \le \min_{x \in X} \{ \mathcal{F}_i(x) \} \\
& \quad + \min_{x \in B_r, w_0 \in \Gamma_0} \{ \mathcal{H}(x, w_0) \} + \theta(\tau_m - \tau_1),
\end{aligned}
\qquad (18)
$$
$$ \min_{x \in B_r, w_0 \in \Gamma_0} \{\mathcal{H}(x, w_0)\} = \min_{w_0 \in \Gamma_0} \{-h^{\eta,g}(x(\tau_1, w_0, \bar{v}))(\tau_m - \tau_1) - g(x(\tau_m, w_0, \bar{v}))\} $$

Then

$$ (\eta + \theta)(\tau_m - \tau_1) + \min_{x \in B_r, w_0 \in \Gamma_0} \{\mathcal{H}(x, w_0)\} + \min_{x \in X} \{\mathcal{F}_i(x)\} \quad (19) $$

is $\varepsilon$-optimal value at $(\tau_1, w_0)$ for $\varepsilon = 2\eta + \theta$ with

$$ g(w) = \int_{\Omega(v_{\min})} L(y, \bar{u}(y), \nabla \bar{u}(y))dy, \quad w \in \Gamma_0, $$

where $\bar{u}$ is a solution to (4) for $\Omega(\bar{v})$.

*Proof.* From the formulae (17), (18) we infer

$$ \min_{x \in X} \{\mathcal{F}_i(x)\} + \min_{x \in B_r, w_0 \in \Gamma_0} \{\mathcal{H}(x, w_0)\} - 2\eta(\tau_m - \tau_1) $$

$$ \leq \max_{x \in B_r} \left\{ \int_{\tau_1}^{\tau_m} \left( - \int_{\Gamma_0} \hat{L}(x(\tau, w, v))dw \right) d\tau - \int_{\Omega(v_{\min})} L(y, \bar{u}(y), \nabla \bar{u}(y))dy \right\} $$

$$ \leq \max_{x \in X} \{\mathcal{F}_i(x)\} + \max_{x \in B_r, w_0 \in \Gamma_0} \{\mathcal{H}(x, w_0)\} + \eta(\tau_m - \tau_1) + \theta(\tau_m - \tau_1). $$

Next, using the definition of value function (5), we get (19).

### 3.3 The algorithm for numerical solution of (P)

In the previous section the last theorem allows us to estimate an $\varepsilon$-optimal value of function (see (6)) for problem (P). As can be seen from the formulas (18) and (19) the essence of the approximation is to be able to calculate the value of the following expressions:

$$ \sup_{x \in X} \left\{ - \sum_{i \in \{1, \dots, m-1\}} x^i (\tau_m - \tau_i) \right\}, \quad \inf_{x \in X} \left\{ - \sum_{i \in \{1, \dots, m-1\}} x^i (\tau_m - \tau_i) \right\}. $$

To achieve this aim we construct a particular directed weighted graph $G$, in which the weight of every edge is the value of the expression $x^i(\tau_m - \tau_i)$. This graph has following properties

1. Every path has length of $m-1$ edges.

2. Every two vertices connected by an edge correspond to points $(\tau, x_1)$ and $(\tau + \Delta\tau, x_2)$ such, that the point $x_2$ is reachable from the point $x_1$ in the next unit of time $\tau + \Delta\tau$.

Therefore by identifying in the graph $G$ the path with lowest (greatest) cost we find the value of the expression $\inf_{x \in X} \{\cdot\} (\sup_{x \in X} \{\cdot\})$.
**The algorithm for generation of the graph G**

1. Let $B$ – a set of trajectories – be a finite set Bezier curves.

2. Let $P$ be a set of points. At the beginning, $P$ contains only one point: $p = w_0$, where $w_0$ is any but fixed point such that $w_0 \in \Omega_0$

3. Create in graph $G$ node which corresponds to point $p$.

4. Calculate $F_g(x_0)$ (equation (8)) where $x_0 = (0, w_0)$, i.e. find a Bezier curve $\alpha \in B$ which minimize value of $G_g(x_0, \cdot)$.

5. For $t = d\tau, \dots, 1$ repeat

(a) Let $P'$ be an empty set of points.

(b) For each point $p$ from $P$ repeat

i. For each Bezier curve $\beta$ from $B$ repeat

A. Find point $p'$ reachable from $p$ under ,,control" $\beta$ in time $t$.

B. Calculate $F_g(x')$ (equation (8)) where $x' = (t, p')$, i.e. find a Bezier curve $\alpha \in B$ which minimize value of $G_g(x', \cdot)$.

C. Create in graph $G$ node which correspond to point $p'$.

D. Create in graph $G$ edge $e_{(p,p')}$ from point $p$ to $p'$.

E. Label edge $e_{(p,p')}$ with a weight $x^t(1-t)$ calculated basing on the indexes $j$ of sets $P_j^{\eta,g}$ which contain points $(t-dt, p)$ and $(t, p')$ generated in the $t$-th step (depending on the difference in those indexes, $x^t$ itself is equal to $-\eta$, 0 or $\eta$).

F. Save $p'$ in $P'$.

(c) Replace set $P$ by $P'$.

6. Generation of the graph $G$ is complete.

**References**

1. Fleming W.H., Rishel R.W., *Deterministic and Stochastic Optimal Control*, New York, Springer-Verlag, 1975

2. P. Fulmański, A. Laurin, J.F. Scheid, J. Sokołowski, *A Level Set Method in Shape and Topology Optimization for Variational Inequalities*, International Journal of Applied Mathematocs and Computer Science, 17, 2007, 413-430.

3. P. Fulmański, A. Nowakowski, J. Pustelnik, *Dynamic programming approach to structural optimization problem – numerical algorithm*, submitted for publication

4. J. Haslinger, R. Mäkinen, *Introduction to Shape Optimization. Theory, Approximation and Computation*, SIAM Publications, Philadelphia, 2003.

5. I. Hlavaček, J. Haslinger, J. Nečas, J. Lovíšek, *Solving of variational Inequalities in Mechanics*, Mir, Moscow, 1996 (in Russian).

6. S. Hüber, G. Stadler, B. Wohlmuth, *A Primal-Dual Active Set Algorithm for Three Dimensional Contact Problems with Coulomb Friction*, SIAM J. Sci. Comput, 30 (2), 2008, 572-596.

7. J. Pustelnik, Approximation of optimal value for Bolza problem (in Polish), Ph.D. thesis, 2009.

8. G. Stadler, *Semismooth Newton and Augmented Lagrangian methods for a Simplified Friction Problem*, SIAM Journal on Optimization 15 (1), 2004, 39-62.