On dual dynamic programming in shape optimization of coupled models

We propose a new method for analysis of shape optimization of coupled models. The framework of the dual dynamic programming is introduced for a solution of the problems. The shape optimization of coupled model is formulated in terms of characteristic functions which define the suport of control. The construction of ε-optimal solution of such a problem can be obtained by solving the sufficient ε-optimality conditions.


Introduction
We consider an optimum design problem with the aim to determine the best location of hollows γ in a given bounded subdomain ω surrounded by the exterior subdomain in R 2 with smooth boundaries ∂ω, = ∂ \∂ω. In the interior subdomain ω the physical phenomenon are described by the linear PDE and in the exterior domain the processes are governed by nonlinear PDE subject to some external function. We took that problem from the paper (Szulc and Zochowski 2015) where more details is presented. The design problem can be considered as a two level optimization problem. At the first level an optimal control within a set of admissible controls is determined for a given location of the source. At the second level an optimal location of the source in terms of its characteristic function is selected in such a way that the resulting value of the cost functional is the best possible within the set of admissible locations. The problem considered in this paper is for an elliptic equation. The problem at the second level is nonconvex, which leads to well known difficulties with the solution procedure. In particular, such difficulties Marta Lipnicka marta@math.uni.lodz.pl Andrzej Nowakowski annowako@math.uni.lodz.pl 1 Faculty of Math, Computer Sciences, University of Lodz, Banacha 22,Lodz,Poland and possible relaxation procedures are discussed e.g. in Kohn and Strang (1986) and Murat and Tartar (1985) mainly from the point of view of existence of solutions in optimum design problems. The standard techniques in classical optimal control theory are based on the lower semicontinuity of some physical quantity (functional) with respect to control and on the compactness of the set of admissible controls. In the optimum design problem the location of the source is optimized. Thus, the lower semicontinuity of the shape functional is required with respect to some family of sets. So, except for very particular cases, there is no optimal location in optimum design problem (see Delfour and Zolesio 2001). That is why many authors confine themself only to numerical analyzis and this is just done in Sokolowski and Zochowski (1999) applying topological derivative and in Szulc and Zochowski (2015). In the last decade some optimum design problems were investigated for time-dependent state equations, including the wave equation or the heat equation (see e.g. Hebrard and Henrot 2003;Münch 2008;Periago 2009;Nowakowski and Sokolowski 2012). Our aim is also to apply the particular variational method, widely used in the classical optimal control theory, which is the dynamic programming technique. In this way, we need neither relaxation nor homogenization of the problem under investigation. We also develop a new computational technic basing on dual dynamic method to characterize the best approximate location of the source explicitly, which is useful for possible applications. In the paper we assume that the number of hollow voids γ is finite but not fixed and non zero i.e. γ = γ 1 ∪...∪γ n , n ∈ N, γ -open ω and such that 2 ≥ vol γ ≥ 1 > 0, with given 1 , 2 . Moreover we assume that all boundaries of γ are smooth. Put D = ∪ ω.
Thus, consider the optimization of the shape and a location of the source γ in ω ⊂ D in the following Problem (P) : where ω D is an open set, such that 2 ≥ vol ω ≥ 1 > 0. The numbers δ, 1 , 2 are given; z d : → R is given target functions, ϕ : ω → R, a fixed function, u : D\γ → R is the state, . Let us observe that by standard elliptic theory for each open set γ ⊂ ω with smooth boundary there exists a unique solution u ∈ H 2 ( ) of (2)-(5). In the problem (P) we are interested not only in a shape of γ but also in a location of γ in ω.
Our main aim is to be able to characterize the optimal (ε -optimal) pair (ū,γ ) and the ε-minimal value of J . Note, thatγ is an optimal (or ε -optimal) open set in our notation. We call a pair (u(·), γ ) "admissible" if it satisfies (2)-(5) and the conditions imposed on γ . Then the corresponding trajectory u(·) is said to be admissible and γ is an admissible set. The set of all admissible pairs is denoted by Ad.
The problems with unknown a subset γ appear in many papers. In Münch (2008) the optimum design for twodimensional wave equation is studied, in Münch (2009) an optimal location of the support of the control for onedimensional wave equation is determined, in Hebrard and Henrot (2005) and Hebrard and Henrot (2003) the optimal geometry for controls in stabilization problem is considered. In all mentioned papers different approaches to the design problems have been investigated, and some numerical results are presented. The existence of optimum design is essential if we have not at hand any sufficient optimality conditions. From the beginning of the last century, under strong influence of Hilbert, the existence issue became one of the fundamental questions in many branches of mathematics, especially in calculus of variations as well as in its branch, the optimal control theory. Of course, following the existence proof, the next step is derivation of necessary optimality conditions and the evaluation of the minimum argument. However, it should be pointed out that for many variational problems the existence of a solution accompanied by some necessary optimality conditions are not sufficient to find the argument of minimum in practice. On the other hand, having at hand a stronger result, i.e. the sufficient optimality conditions for a minimum in a specific problem, replaces the requirement for the existence.
In the present paper the framework of dynamic programming together with sufficient optimality conditions (the so-called verification theorem for relative minimum) is proposed for a solution of the optimum design problem. Different approaches are given in Sokolowski and Zochowski (1999), Bednarczuk et al. (2000), Nowakowski and Sokolowski (2012), Fulmański et al. (2014), and Fulmański and Nowakowski (2014). In Sokolowski and Zochowski (1999), Bednarczuk et al. (2000), and Nowakowski and Sokolowski (2012) the shape problem is formulated in terms of characteristic functions which define the support of control while in Fulmański et al. (2014) and Fulmański and Nowakowski (2014) also the framework of dynamic programming together with sufficient optimality conditions is proposed for a solution of the optimum design problem. In Fulmański et al. (2014) the problem of dividing tube is investigated by reduction of the shape optimization problem to a classical control problem. To do that the authors apply the level set approach to build deformation, following Zolésio (Sokolowski and Zolesio 1992) they introduce a field depending on control, which define the type of deformation, and formulate the control problem governed by ordinary differential equation (defined by field). Next the classical dynamic programming is developed for such a problem. In Fulmański and Nowakowski (2014) a little similar approach by level set method is applied to the optimization problem of deformable structures described by Navier-Stock equation but then the dual dynamic programming is delveloped to derive the verification theorem (again for dynamics governed by ordinary differential equation). In both papers completely different numerical algorithms are constructed the base of which are different approximate verification theorems. However we should stress that approaches based on the level set method described in Fulmański et al. (2014) and Fulmański and Nowakowski (2014) use smaller set of admissible deformation (described by a field depending on control) i.e. in our case a significantly smaller set of the sources γ . In the presentation (Kaźmierczak 2008) similar ideas to those of Fulmański et al. (2014) and Fulmański and Nowakowski (2014) are mentioned but nothing is precisely formulated and proved. Our goal is not the standard analysis of the problem as e.g. in Delfour and Zolesio (2001) or Fulmański et al. (2014) and Fulmański and Nowakowski (2014) but the approximate solution by application of the sufficient ε -optimality conditions given by the dual dynamic programming directly to our problem (P). That means we admit full set of admissible sources γ . This approach seems to be new and the result obtained is original, to our best knowledge.
We provide a dual dynamic programming approach to control problems (2)-(5). This approach allows us to obtain the sufficient conditions for optimality in problem considered. We believe that the conditions for problems of type (P) in terms of the dual dynamic programming, that we formulate here, have not been provided earlier.
There are two main difficulties that must be overcome in such problems as (P). The first one consists in the following observation. We have no possibility to perform perturbations of the problem -as it is considered in the fixed set D with boundary condition (5) -which can be compared to the one-dimensional case given in Bellman (1957) and Fleming and Rishel (1975). The second one is that we deal with elliptic equation for state and controls (2 ). The technique we apply is similar to the methods from Galewska and Nowakowski (2006) and Nowakowski (1992). The main idea of the methods from Galewska and Nowakowski (2006) and Nowakowski (1992) is that they carry over all objects used in dynamic programming to dual spacespace of multipliers (similar to those which appear in the Pontryagin maximum principle). Next, instead of classical value function (which for problem (P) makes no sense), we define an auxiliary function V (x, p) satisfying the second order partial differential equation of dual dynamic programming (compare Galewska and Nowakowski 2006). Investigations of the properties of this function lead to an appropriate verification theorem. We introduce also the concept of an optimal dual feedback control and provide new sufficient ε -optimality conditions determined within our framework. Just by using the approximate differential equation of dual dynamic programming (9 ) (see below) and ε-optimal dual feedback control (section 4) we are able to solve problem (P) completely, from the approximate point of view i.e we find aproximate location of γ ε as well for ε-optimal value of J .

Dual dynamic programming approach
In Galewska and Nowakowski (2006) a new approach to dynamic programming was developed using some ideas of Nowakowski (1992). In Nowakowski (1992), instead of considering notions of dynamic programming such as value function S(t, s) or Hamilton-Jacobi equation in the space (t, s), a new space -the dual space is proposed and new notions of dual dynamic programming are defined: an auxiliary function, a dual optimal value and a dual Hamilton-Jacobi equation which the auxiliary function should satisfy. The dual space in Nowakowski (1992) is, in fact, defined by conjugate (dual) functions (variables) which appear in Pontryagin maximum principle. In turns out that this approach works also in control problems governed by elliptic equations (see Galewska and Nowakowski 2006). That means: in dual approach to dynamic programming the perturbation of optimal value is not needed -instead we deal with an auxiliary function. However, there is a price to be paid for that, as we have to impose on the auxiliary function some additional condition, called the transversality condition. Our aim is to apply ideas from Galewska and Nowakowski (2006) to the optimum design problem (P). To this end we need to define dual notions in some dual space. Thus let P ⊂ R 2+2 be an open (dual) set of the variables (x, p) = (x, y 0 , y), x ∈ D , y 0 < 0, |y| ≥ 0. We shall also use the sets Denote by W 2:3 (P ) the specific Sobolev space of functions of two variables (x, p) having up to the second order weak or generalized derivatives (in the sense of distributions) with respect x and up to the third order weak derivatives with respect to the variable p. Our notation for the function space is used for the functions depending on the primal variable x, and the dual variable p, the primal and dual variables are independent and the functions in the space W 2:3 (P ) enjoy different properties with respect to x and p. Let V (x, p) of W 2:3 (P ) be a (auxiliary) function defined on P and satisfying the following condition: for (x, p) ∈ P . Here, V y 0 , V y , and V p denote the partial derivatives and the gradient with respect to the dual variables y 0 , y, and p = (y 0 , y), respectively. Now, we denote by p(x), x ∈ D\γ , the dual trajectory, while u(x), x ∈ D\γ stands for the primal trajectory. Let us put Using the function u it is possible to come back from the dual trajectories p(x), x ∈ D\γ, lying in P to the primal functions u(x), x ∈ D\γ . How to find V y is precisely described below. Further, we confine ourselves only to those admissible trajectories u(·), for which there exist functions Actually, it means that we are going to study problem (P) possibly in some smaller set Ad u , which is determined by the function (7). Next, we define a dual optimal value S u D for the problem (P) by the formula In order to prove the verification theorem we require that the function V (x, p) satisfies the second order partial differential equation in dynamic programming form: with boundary conditions Let us note that in (9) maximum is taken instead of supremum since it is attained. We shall not discuss here the question of existence of solution to (9) and satisfying condition (6). We simply assume in verification theorem (in next section) that such a function exists.

A verification theorem
In this section we formulate and prove our main theorem called "verification theorem", which provides the sufficient optimality conditions for (P) in terms of a solution V (x, p) of the second order partial differential equation of dynamic programming (9). Let us fixȳ 0 < 0 . Define the set Theorem 1 Assume that there exists a W 2:3 (P ) solution V (x, p) of (9) on P with boundary condition (10) such that (6) holds and letū(x, Then (u(·),γ ) is an optimal pair relative to all (u(·), γ ) ∈ Adū.
Proof Our proof starts with the observation that from transversality condition (6), we see that for x ∈ D\γ , Now define a function D\γ x → W (x, p(x)) by We conclude from (12)- (14) that Hence, (9) and (15) imply and finally, after integrating (16) and applying (14), we have Thus from (17) and the Green formula it follows that where ν(·) is the exterior unit normal vector to ∂ . In the same manner applying (11) and (14) we have Combining (18) with (19) gives which completes the proof.

Optimal dual feedback control
This section is devoted to the notion of an optimal dual feedback control. We present appropriate definitions and the sufficient conditions for optimality which follow from the verification theorem.
In addition, the optimal value S u D (see (8)) is defined by the pair (ū(·),γ ): and there is V (x, p), (x, p) ∈ P belonging to W 2:3 (P ), satisfying (6), such that x V y 0 (·, p(·)) ∈ L 2 ( ) and y 0 x V y 0 (s, p(s))ds = −S u D , V y (x, p) = −u(x, p). Now, we formulate and prove the sufficient optimality conditions for the existence of an optimal dual feedback control, again in terms of the auxiliary function V (x, p).

Approximate optimality for the problem (P)
In the last two subsections we discussed optimality conditions in terms of dual dynamic programming. In practice it is difficult (or even impossible) to solve equations stated there in exact form. In fact we solve such a system using different approximate (numerical) methods. Therefore what we can get then that is eventually approximate optimality. This is why in this section we present dual dynamic approach to sufficient conditions for approximate (ε-optimality) optimality. Just the dual ε-optimality conditions are the base for construction of computational, approximate optimality. Let us recall that for given u and y 0 (see (8)) optimal value is defined as Dual ε-optimal value for problem (P) we call each value Let us fix m > 0 and put P = (x, y 0 , y); y ∈ ⊂ R 2 , y 0 < 0, y > 0, As for ε-optimal value we use inequality it suggests that expressions allowing to derive Theorem 1 should satisfy also suitable inequalities. Thus we shall use the following inequalities for auxiliary functionṼ : the dual Hamilton-Jacobi inequality with boundary conditions Since we want to apply our theory to numerical solutions of (2)-(5), therefore instead of the equation we shall deal with the inequality: Thus in this section by the set of admissible sets and states i.e. satisfying (28) we denote Ad ε .

Scalability and convergence of the algorithm
First we would like to stress that using any package to solve nonlinear differential equation we are able only to find approximate solution i.e. the numerical solutions do not satisfy in our case the Eqs. 2, 9. Therefore any sufficient optimality conditions for approximate optimality should take into account those circumstances. In our algorithm we do it, i.e. we solve with any computational package the inequalities (28), (26) and check whether the calculated solutions satisfy (28), (26) with given ε. The number N in points 2., 3. of the above Algorithm may be arbitrarily large. We can then divide N on smaller parts and do calculations for (28) and point 3. (on those parts) independly (on different processors) and next to find minimal value of J as a minimum of minimal values of J on those parts. If we were looking for γ being balls (or systems of balls) then we can cover the whole bounded domain ω with finite number of balls of any positive fixed radius. Then the convergence of the algorithm is obvious. However if we are interested not only in locations but also in shapes of γ then from the theoretical point of view the convergence of algorithm is still almost obvious as the closure of ω is a compact set, thus there is always finite covering of the closure of ω by open sets γ . Hence we can calculate at least theoreticaly the best of them using te above algorithm and any other set γ , by continuity of our functional is near one of a set from the covering family.

The parameters
In the numerical simulations we take Fig. 1).
-three balls where (x 1 1 ,x 1 2 ) -the center of γ 1 , where (x 2 1 ,x 2 2 ) -the center of γ 2 , We assume that ω ε = ω \ γ , ε = 0.1,γ -the optimal γ , y 0 ε = −0.005 and N = 500. For a randomly generated control, we calculate the value of J . For this purpose, we define the function z d (or z d1 and z d2 respectively) by defining the domain for each type of control. It means that we define We are looking for optimal control as ball (or systems of balls) that generate minimal values of J and satisfy (34)-(35).

Numerical calculations
In order to make calculations we apply Matlab2016a and implement in this application the following steps: 1. Let ε = 0.1. 2. Let fixed number of balls (one, two or three). 3. We calculate function z d (or z d1 and z d2 respectively). z d is a solution of the following boundary value problem: where f zd (x) = x 1 for x ∈ D\γ zd when we consider one ball. 4. If we consider two or three balls we solve above boundary value problem due to two functions f zd1 (x) = x 1 + x 2 for x ∈ D\γ zd12 and f zd2 (x) = x 1 − x 2 for x ∈ D\γ zd12 . 5. For given parameters we randomly generate γ for fixed ε. 6. For given γ , ω and we solve the following systems of equations where ω .

Numerical examples
We present four examples, each for an otherwise defined control.

Example 1
Let us γ as one ball. For given parameters we realize numerical calculations in the following way: -We fixed ε = 0.1.
-For given parameters we repeat 3. -5. N = 500-times. We select the minimum value among all values from (1). It means that we find minimal value of J (γ n ), n = 1, ..., N and corresponding to it pair denote by -We generate the set Y -We calculate the functionṼ from (26)-(27).

Example 2
Let us γ as two holes defined as two balls. For given parameters we realize numerical calculations in the same way like in Example 1. We obtain the following result: -Functions z d1 and z d2 (Figs. 6 and 7).

Example 3
Let us γ as three holes defined as three balls. For given parameters we realize numerical calculations in the same way like in Example 1. We obtain the following result: -Functions z d1 and z d2 (Figs. 9 and 10).  Let us γ as one hole defined as non-convex set -two balls with nonempty intersection.
Remark 2 In Example 4 we considered the controls as nonconvex holes. To calculate function z d we defined domain with one hole defined as a ball. For this case the verification theorem is satisfied.

Remark 3
The approximate value of our functional J is similar like the result in the numerical example presented in Szulc and Zochowski (2015). However we know how far we are from the inf J (γ ), while in Szulc and Zochowski (2015) is only known some step of calculation.

Conclusions
The dual dynamic programming by construction furnishes the sufficient optimality conditions, therefore it is a powerful tool for solution of optimization problems which enjoy the special structure. This is especially important when we deal with approximate (numerical) solutions. In the paper for the optimum design problem the dual dynamic programming and dual feedback notions are developed as well their approximate counterparts. Next the solutions of that problem are characterized in terms of verification theorems. The approximate verification theorem allows us to calculate approximate solution but also to state how far we are from the inf J (γ ). In Szulc and Zochowski (2015) the value of calculated numerically functional is almost the same (we took the same parameters in our example) but the authors in Szulc and Zochowski (2015) can only assert that it is some step of calculation of the minimal value. In a subsequent paper we are going to apply the dynamic programming technique to more general shape optimization problems.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.