∎

11institutetext: C. Guerrero (Corresponding author), I. Lera, C. Juiz
22institutetext: Computer Science Department, University of Balearic Islands,
Crta. Valldemossa km 7.5, E07122 Palma, Spain
Tel.: +34-971172965
Fax: +34-971173003
C. Guerrero
22email<EMAIL_ADDRESS>
I. Lera
22email<EMAIL_ADDRESS>
C. Juiz
22email<EMAIL_ADDRESS>.

# Genetic algorithm for multi-objective optimization of container allocation
in cloud architectures

Carlos Guerrero  Isaac Lera  Carlos Juiz

(Received: date / Accepted: date)

###### Abstract

The use of containers in cloud architectures has become widespread because of
advantages such as limited overhead, easier and faster deployment and higher
portability. Moreover, they are a suitable architectural solution for
deployment of applications created using a microservices development pattern.
Despite the large number of solutions and implementations, open issues have
not been addressed in container automation and management. Container resource
allocation influences system performance and resource consumption so it is a
key factor for cloud providers. We propose a genetic algorithm approach, using
the Non-dominated Sorting Genetic Algorithm-II (NSGA-II), to optimize
container allocation and elasticity management due to the good results
obtained with this algorithm in other resource management optimization
problems in cloud architectures. The optimization has been focused on a tight
use of the resources and a reduction of the network overhead and system
failure rate. A model for cloud cluster, containers, microservices and four
optimization objectives is presented. Experimental results have shown that our
approach is a suitable solution to address the problem of container allocation
and elasticity and it obtains better objectives values than the container
management policies implemented in Kubernetes.

###### Keywords:

Cloud containers Microservices Resource allocation Genetic algorithm Multi-
objective optimization Performance evaluation

##  1 Introduction

Microservices architecture is a development pattern that has grown in
popularity in recent years. This pattern defines an application as a set of
independent small and modular services executing a single task and the
specification of the interoperability of these microservices to achieve the
application requirements. Several companies have adopted this practice and
traditional service cloud applications are being migrated to microservices
balalaie16microservices  for all their benefits: making easier to ship and
update applications, allowing independently update and redeploy of parts of
the application, getting closer development and operations teams, allowing
continuous release cycle, simplifying orchestration of applications across
heterogeneous cloud data centers and so on fazio16open .

The containerization of applications is one of the enabling technologies for
microservices architectures. Containerization is a method for operating system
(OS) virtualization. Virtual instances of OS, called containers, are a form of
isolation for the OS environment and its file system, and they are executed
over a single host and a single kernel. Therefore, each microservice of an
application is deployed in a container, without launching a complete virtual
machine. By this, system overhead is reduced as security and isolation of the
environments are maintained. The deployment of a microservice is as simple as
starting the execution of a new container of the microservice. Therefore,
microservices can be scaled up by just creating new containers until the
desired scalability level is achieved. Docker is one of the most successful
implementations of container architectures peinl2016docker . Container
clusters have been developed to easily deploy containers in the cloud and to
enable features as consolidation, elasticity, load balancing and scalability
in container architectures. Docker Swarm, Apache Mesos and Google Kubernetes
are examples that provide automatic deployment and orchestration
peinl2016docker .

Despite the fast growth of these technologies and the existence of de-facto
solutions, there are important open issues in scheduling and management of
microservices in the cloud. Allocation of containers across the physical nodes
of a cluster influences system performance and reliability, as well as
elasticity, and its automatic management, influences application performance.
The problem of elasticity in virtual machines environments has been widely
studied caballer2015dynamic  but, to the best of our knowledge, not in the
case of container-based clouds. Current container cluster managers, as
Kubernetes, implement simple approaches for the automatic scalability of
containers and the allocation of containers to physical machines, only focused
on physical resource usages and general threshold values fazio16open . We
propose to improve the elasticity schedulers and runtime container allocators
to obtain better results in terms of system provisioning, system performance,
reliability and network overhead, four of the traditional research challenges
in cloud architectures zhang2010cloud .

Resource management optimization is a NP-complete problem and it needs to be
addressed by metaheuristic approaches wei2010game . The optimal solution can
be only obtained with an evaluation of all the possible combinations. The
literature reflects that evolutionary approaches, as Genetic Algorithms (GA),
are common solutions to resource management in cloud zhan2015cloud .
Furthermore, Non-dominated Sorting Genetic Algorithm-II (NSGA-II) is one of
the most usual solutions for multi-objective optimization vonlucken2014survey
.

We propose to use NSGA-II to impolement the container allocation strategy and
the automatic elasticity management (scalability level of the microservices)
by the optimization of four objectives:

  * •

the provisioning of new applications and the elasticity of the currently
deployed ones by keeping a uniform distribution of the workload across the
cluster;

  * •

the performance of the deployed applications and their assigned resources by
considering a suitable scalability level of their microservices and a uniform
distribution of the microservices workload across their containers;

  * •

the reliability of the microservices by considering a suitable scalability
level and a correct distribution avoiding single points of failure;

  * •

the network overhead of the communications between microservices by placing
the containers of related microservices in physical machines with short
network distances.

We summarize the main contributions of our work in:

  * •

The improvement of resource management in container clusters by optimizing
container allocation to physical machines and the automatic container
scalability;

  * •

The use of four decision criteria —microservices reliability, microservices
network overhead, balanced use of the physical machines and balanced workload
over the containers— for the resource management optimization;

  * •

The implementation and the study of an evolutionary approach, the genetic
algorithm NSGA-II, to solve this type of multi-objective optimization problem;

  * •

The validation of our approach with workload and system models based on data
from real deployed scenarios .

##  2 Related work

Resource management optimization is a common research topic in the field of
Cloud Computing manvi2014resource . Several surveys have been published
including a wide range of techniques to address resource management problems
as, for example, scheduling singh2016survey , provisioning sukhpal2016cloud
and allocation beloglazov2012energy ; akhter2016energy .

In addition, evolutionary approaches and genetic algorithms are commonly used
in resource management in cloud environments. The adequacy of genetic
algorithms in general, and NSGA-II in particular, has been widely demonstrated
along the research bibliography zhan2015cloud . Most of the researches are
addressed to improve the placement of virtual machines in physical machines
based on objectives from very different points of view, e.g. performance,
reliability or energy efficiency. For example, Pascual et al. studied three
multi-objective optimization algorithms, including NSGA-II, and demonstrated
that the use of evolutionary approaches is a low-cost optimization solution to
obtain better performance and to reduce the energy consumption
pascual2015towards .

There is a wide range of implementations of container-based solutions but
issues as container orchestration and self-management remain open
peinl2016docker . The number of studies for resource management in container-
based architectures is very limited. To the best of our knowledge, our
solution is the first approach that uses genetic algorithms and only three
recent works deals with container allocation problem: two of them based on
Linear Programming Solvers and another one based on game theoretic problem.

Guan et al. guan2016application  proposed an application-oriented algorithm
that decides the physical node where a new container is scheduled by
minimizing the application deployment cost. Kang et al. also implemented a
container scheduler for new arrivals of containers but considering the
objectives of saving energy consumption and guaranteeing the performance level
kang2016cloud . These two first works differs from our proposal because they
solve the allocation problem only for new container arrivals while the current
allocated containers are not migrated between machines.

Finally, Xu et al. presented an approach based on the Stable Marriage Problem,
a many-to-one matching problem, to reduce the response time of costumers’ jobs
and to improve the physical computational resources utilization xu2014novel .
In this approach, not only the optimization objectives, , or satisfaction as
they called it, need to be measured but also a rank of physical machines for
each container and a rank of containers for each physical machine need to be
calculated. These rank functions are based on criteria as selecting the
machine with more resources, or the proportionality between container resource
consumption and machine capacity. When a physical machine and a container are
the best ranked between themselves, the container is allocated in the machine.
The satisfaction function is used after each new matching to ensure that the
allocation improves the optimization objectives. The allocation of the
machines is driven by the rank functions instead of the optimization
objectives. In our approach, rank functions are not considered, and our
objective functions can not be directly translated to rank functions.

Our proposal is to use NSGA-II, a general accepted approach for multi-
objective optimization problems in cloud architectures, in a resource
management problem, container allocation and elasticity that, from our current
knowledge, has not been previously addressed with evolutionary strategies.

##  3 System model

Our approach needs to use a formal definition of the system by modeling
applications, microservices, physical machines and network. For a better
understanding, Table 1 summarizes the parameters of the system model described
in the remainder of this section.

We consider a set of applications $A$ following a development pattern based on
microservices. Each of these applications ($app_{j}$) is characterized with
the number of user requests ($ureq_{j}$) and the microservices stack. The
microsevices stack is the set of microservices and the interoperability
between them that implements the requirements of an application. These
interoperability relationships are established when microservices consume
results generated by other microservices. By this, the microservices stack is
modeled as a directed graph. The nodes are the microservices ($ms_{i}$) and
the edges represent the interoperability between them. Two nodes are connected
by an edge, $(ms_{provider},ms_{consumer})_{prov/cons}$, if the target
microservice $ms_{consumer}$ consumes the results of the origin one
$ms_{provider}$.

Each microservice is characterized as a tuple
$(msreq_{i},res_{i},thr_{i},fail_{i})$ where, $msreq_{i}$ is the number of
microservice requests needed to satisfied one application user request
($ureq_{j}$); $res_{i}$ is the computational resources consumed to satisfied
one request of the microservice; $thr_{i}$ is the threshold level for the
resources consumption above which the service performance is degraded and the
microservice results on a bottleneck for the application $app_{j}$; and
$fail_{i}$ is the failure rate of the microservice. The values for $ureq_{j}$
and $msreq_{i}$ depend on the workload of the system and the values for
$res_{i}$, $thr_{i}$ and $fail_{i}$ depend on the implementation of the
microservices.

Each microservice is executed in the system encapsulated in one or more
containers $cont_{k}$, represented as $ms_{i}\equiv cont_{k}$. The number of
containers depends on the scale level of the microservice ($scale_{i}$). The
microservice resource consumption is divided uniformly among the containers,
so the computational resource consumption of the container $res_{k}$ is
calculated as $\frac{ureq_{j}\times msreq_{i}\times res_{i}}{scale_{i}}$.

The set of containers ($C$) is executed in a container cluster with a set of
physical machines $pm_{l}$. The relationship $alloc(cont_{k})=pm_{l}$
indicates that container $cont_{k}$ is allocated to the physical machine
$pm_{l}$. The notation can be summarized to represent the physical machine
that allocates a container of a microservice $ms_{i}$ using
$alloc(ms_{i})=pm_{l}$ if $cont_{k}\equiv ms_{i}$ and
$alloc(cont_{k})=pm_{l}$.

Each physical machine is characterized by the tuple $(cap_{l},fail_{l})$ where
$cap_{l}$ is the computational capacity and $fail_{l}$ is the failure rate of
the node. A constraint of the system is that the sum of computational
resources used by the containers allocated to a physical machine needs to be
lower than the computational capacity of this physical node. Finally, the
physical machines are interconnected with a physical network $N$, where the
paths between nodes $pm_{l},pm_{l^{\prime}}$ are characterized by their
network distance $dist_{pm_{l},pm_{l^{\prime}}}$.

Table 1: Summary of the system model Element | Parameter | Description
---|---|---
Application | $app_{j}\in A$ | application with id. $j$
| $|A|$ | total number of applications in the cluster
| $ureq_{j}$ | number of user request for application $app_{j}$
Microservice | $ms_{i}\ inS$ | microservice with id. $i$
| $|S|$ | total number of microservices
| $(ms_{i^{\prime}},ms_{i})_{prod/cons}$ | microservice $ms_{i}$ consumes microservice $ms_{i^{\prime}}$
| $msreq_{i}$ | microservice requests number $ms_{i}$
|  | needed for each $ureq_{j}$ request
|  | from application $app_{j}$
| $res_{i}$ | computational resources required for a
|  | microservice request
| $thr_{i}$ | threshold for normal operation of a
|  | microservice $ms_{i}$
| $fail_{i}$ | failure rate for a microservice $ms_{i}$
Container | $cont_{k}\in C$ | container with id. $k$
| $|C|$ | total number of containers
| $cont_{k}\equiv ms_{i}$ | container $cont_{k}$ encapsulate/execute
|  | microservice $ms_{i}$
| $scale_{i}$ | number of containers for $ms_{i}$
| $res_{k}$ | computational resource consumption of
|  | container $cont_{k}$
Physical machine | $pm_{l}\in P$ | physical machine with id. $l$
| $|P|$ | total number of physical machines
| $alloc(ms_{i}/cont_{k})=pm_{l}$ | physical machine $pm_{l}$ allocates
|  | service $ms_{i}$ / container $cont_{k}$
| $cap_{l}$ | computational capacity of physical
|  | machine $pm_{l}$
| $fail_{l}$ | failure rate for a physical
|  | machine $pm_{l}$
Network | $dist_{pm_{l},pm_{l^{\prime}}}$ | network distance between $pm_{l}$ and $pm_{l2^{\prime}}$

###  3.1 Optimization objectives

In our approach, the optimization of the objectives is achieved through the
allocation of the containers and the definition of the scale level, i.e.,
relationship $allocation(cont_{k})$ and $scale_{i}$ are determined $\forall
cont_{k}$ and $\forall ms_{i}$ using a set of optimization objectives. The
objectives to optimize are: (OBJ-i) the container workload avoiding the
presence of bottlenecks; (OBJ-ii) the balanced used of the cluster to
facilitate the future admission and provision of new applications; (OBJ-iii)
the reliability of the applications by an evenly distribution of the
containers through the cluster nodes; and (OBJ-iv) the intercommunication
overhead by allocating related microservices to physical machines with short
network distances.

For the first optimization objective (OBJ-i), we consider that the workload of
the containers is appropriated when the resource consumption of the container
$res_{k}$ is similar to the resource consumption threshold of the microservice
$thr_{i}$. Thus, a high scalability level of a microservice is penalized when
the container is underused –the container resource consumption is lower than
the microservice resource threshold– and a low scalability level is penalized
when the number of microservice requests is very high –the resource
consumption is higher than the microservice resource threshold–. We define a
metric call Threshold Distance that is the difference between the resource
consumption of a container and the threshold value of a microservice. This is
formalized in Equation 1.

$Threshold\ Distance\ =\sum_{\forall ms_{i}}\left|\frac{ureq_{j}\times
msreq_{i}\times res_{i}}{scale_{i}}-thr_{i}\right|$ (1)

For the second optimization objective (OBJ-ii), we consider that the cluster
is balanced when the use of the physical nodes is uniform across the cluster
—the computational resources consumption is uniformly distributed among the
physical machines—. If a physical machine does not have any container
allocated, it is not considered in this calculation, because it could be
switched off. We propose to use the standard deviation of the percentage of
resource usages of the physical nodes to evaluate the balance of the cluster
xu2016location . This is formalized in Equation 2.

$\displaystyle Cluster\ Balanced\ Used\ =\sigma(PM_{usage}^{pm_{l}},\quad
if\quad\exists\ ms_{i}\quad|\ alloc(ms_{i})=pm_{l})$ (2)

where:

$\displaystyle PM_{usage}^{pm_{l}}=\frac{\sum_{ms_{i}}\frac{ureq_{j}\times
msreq_{i}\times res_{i}}{scale_{i}}}{cap_{l}}\quad\forall ms_{i}\quad|\
alloc(ms_{i})=pm_{l}$ (3)

For the third optimization objective (OBJ-iii), we measure the reliability of
the system using the failure rate of the applications. An application fails
when any of their microservices fails. A microservice fails when all the
container replicas fail. A container fail is generated by a fail in the
container, $fail_{i}$, or by a fail in the physical machine that allocates the
container, $fail_{l}$. This is mathematically modeled in Equation 4.

$\displaystyle System\ Failure\ =\sum_{\forall ms_{i}}Service\
Failure(ms_{i})$ (4)

where:

$\displaystyle Service\ Failure(ms_{i})\ =$ (5) $\displaystyle=\prod_{\forall
pm_{l}\ |\ allocation(ms_{i})=pm_{l}}\left(fail_{l}+\prod_{\forall ms_{i}\ |\
allocation(ms_{i})=pm_{l}}fail_{i}\right)$ (6)

For the fourth optimization objective (OBJ-iv), we approximate the network
overhead between the microservices of an application with the average network
distance between all the pairs of consumer and provider containers.
Consequently, we consider the mean value of the distances between the replicas
of a microservice and the containers executing microservices consumed by the
former. Equation 7 formalizes this metric.

$\displaystyle Total\ Network\ Distance=\sum_{\forall ms_{i}}Service\ Mean\
Distance(ms_{i})$ (7)

where:

$\displaystyle Service\ Mean\ Distance(ms_{i})=$ (8)
$\displaystyle=\frac{\sum_{\forall cont_{k}|cont_{k}\equiv
ms_{i}}\left(\sum_{\forall cont_{k^{\prime}}\equiv ms_{i^{\prime}}\
|(ms_{i^{\prime}},ms_{i})_{prov/cons}}dist_{alloc(cont_{k}),alloc(cont_{k^{\prime}})}\right)}{|cont_{k}|\times|cont_{k^{\prime}}|}$
(9)

To sum up, our aim can be defined as a four-objectives optimization problem
where the solution establishes (a) the scalability level of each microservice
and (b) the allocation of containers to the physical machines, by minimizing
(OBJ-i) Threshold Distance, (OBJ-ii) Cluster Balanced Used, (OBJ-iii) System
Failure and (OBJ-iv) Total Network Distance functions. The problem can
therefore be formulated as:

Determine:

$\displaystyle scale_{i}=|\\{cont_{k}\\}|\quad|\quad cont_{k}\equiv
ms_{i}\quad\forall ms_{i}$ (10) $\displaystyle alloc(ms_{i})\quad\forall
ms_{i}$ (11)

By minimizing:

$\displaystyle\sum_{\forall ms_{i}}\left|\frac{msreq_{i}\times
res_{i}}{scale_{i}}-thr_{i}\right|$ (12) $\displaystyle\sigma(PM_{usage})$
(13) $\displaystyle\sum_{\forall ms_{i}}Service\ Failure(ms_{i})$ (14)
$\displaystyle\sum_{\forall ms_{i}}Service\ Mean\ Distance(ms_{i})$ (15)

Subject to the constraint:

$\displaystyle\sum_{\forall
cont_{k}|alloc(cont_{k})=pm_{l}}res_{k}<cap_{l}\quad\forall\ pm_{l}$ (16)

This optimization problem is a NP-complete problem because all the possible
scalability levels and allocations should be evaluated to find the optimal
solution. A metaheuristic method needs to be used to address the problem
wei2010game .

##  4 Genetic algorithm proposal

Genetic algorithms (GA) are metaheuristic approaches for solving multi-
objective optimization problems. Our solution is based on the Non-dominated
Sorting Genetic Algorithm-II (NSGA-II) that is considered a standard approach
for multi-objective problems. In the field of GA, and along this manuscript,
the following terms are used interchangeably: solution, chromosome and
individual; set of solutions and population; algorithm iteration and
generation; offspring, next population and next algorithm iteration.

GA implementation requires to define the representation of the chromosome, the
crossover and mutation operators, the fitness function, the selection
operator, the offspring generation and the execution parametrization. The
quality of the solutions is influenced by how these elements are defined
mitchell1998introduction . The next subsections explain the details of the
adaptation of each part of the GA to our optimization problem. Additionally,
we include a general structure of the genetic algorithm in Algorithm 1 for a
better understanding.

###  4.1 Chromosome representation

In GA, the solutions for the problem are usually represented using a string-
based notation. Each of these solutions is called chromosome. We define a
string-based notation to represent the number of containers (scalability
level) for each microservice, and the allocation of these containers to the
physical machines (Figure 1).

The representation of the allocation is an array, called microservices array,
where the index is the microservice identifier and the content is a list,
called allocation list, with the identifiers of the physical machines where
the containers for the given microservice are allocated. The number of
elements in the allocation list is the scalability level for a microservice.
The same physical machine identifier can be repeated several times in the
allocation list because the physical machine can allocate several containers
of the same microservice.

Figure 1: Chromosome representation example for a system with 3 physical
machines allocating 6 microservices deployed in 13 containers

###  4.2 Crossover and mutation operators

Each iteration of a GA generates a new population for the given problem. Each
of these new generations tends to approach to better values for the
objectives. This trend is based on the biological evolution concepts of
crossover and mutation, where new individuals are created by mating the best
ones of the previous generations. In GA, crossover and mutation operators
define the way that parent solutions are combined and modified to generate new
ones.

Crossover operator combines two solutions to generate, usually, two new ones.
This combination is defined over the chromosome representation taking
alternatively pieces from the first and the second chromosome. In our case,
the crossover is applied independently to each single allocation list. In each
of those lists, a single-point crossover operator is applied
mitchell1998introduction . This operator combines two chromosome pieces, one
from each parent, to generate the first individual. The second one is
generated by the combination of the opposite pieces of the chromosomes.

The single-point crossover operator generates a random number, $r$, between 1
and the number of elements in the smaller parent allocation list,
$min($$length(c1),$ $length(c2))$. The outputs of the operator are two new
allocation lists obtained from the concatenation of the $[1,r]$ elements of
the first list and the $[r+1,length(c2)]$ elements of the second list —for the
first output— and the concatenation of the $[1,r]$ elements of the second list
and the $[r+1,length(c1)]$ elements of the first list —for the second output—.
A different random value $r$ is generated for each allocation list. Figure 2
shows an example of the crossover operator applied to two parents.

Figure 2: Example of the crossover operator applied to two parent solutions to
obtain two new individuals

Mutation operator is usually applied to just one individual. The objective of
mutation operator is to modify solutions randomly to avoid local solutions in
the solution search space mitchell1998introduction  .

Figure 3: Examples for the three mutation operators

We have defined three mutation operators: Swap Mutation, which shuffles the
array positions (allocation lists) of the microservice array; Shrink Mutation,
which removes a random number of physical machines from all the allocation
lists, i.e., the scalability level is decreased for all the microservices; and
Growth Mutation, which increases the allocation lists by adding a random
number of physical machine identifiers, i.e., the scalability level is
increased for all the microservices. Figure 3 represents examples of the three
mutation operators.

###  4.3 Fitness function and selection operator

The NSGA-II algorithm includes a specific definition of how the fitness
function and the selection operator have to be used. NSGA-II is designed for
multi-objective optimizations. Thus, the solution is not evaluated with a
single fitness value. The suitability of the solution is expressed by a tuple
with the values for each of the optimization objectives. In our case, the
fitness value is represented with the values of the four functions to minimize
(Equations 1, 2, 4 and 7). If the constraints for a solution are not
accomplished (Equation 16), the fitness value is fixed to infinite for all the
criteria.

The selection operator establishes how the parents of a crossover operation
are selected from the population. To explain the details of the selection
operator defined in the NSGA-II, we need to previously introduce three
concepts deb2002fast . The first one is the dominance. A solution $s_{1}$
dominates another solution $s_{2}$ if the fitness values of $s_{1}$ are
better, for all the objectives, to the values of $s_{2}$. In the same way, a
solution $s^{\prime}_{1}$ non-dominates another solution $s^{\prime}_{2}$ if
the fitness values of $s^{\prime}_{1}$ are better for a subset of objectives
and worse to the others. Secondly, the Pareto optimal front is the set of
solutions that are non-dominated by other solutions. Therefore, each solution
in the Pareto front optimizes one or more objectives, but not all of them, in
comparison to the other solutions in the front. Finally, the crowding distance
is a metric to estimate the density of solutions surrounding a given one. The
crowding distance is calculated obtaining the average distance along each of
the objectives of the closest neighbors of the solution.

The NSGA-II algorithm selection operator is based on the idea that the
probability of selecting a solution would be higher when it is a non-dominated
solution and it has less solutions surrounding it. Therefore, the NSGA-II
algorithm sorts the solutions in front levels, and all the solutions in the
same front level are ordered using the crowding distance (lines 19 to 22 in
Algorithm 1). Once all the solutions are sorted, it applies a binary
tournament selection operator over the sorted elements (lines 12 to 13): two
solutions are selected randomly and the first one in the ordered list is
finally selected konak2006multi .

The front levels are generated iteratively. In a first step, the Pareto
optimal front is calculated, and it is considered the first front level. In a
second iteration, all the solutions in this first front level are removed, and
the new Pareto optimal front is calculated, obtaining the second front level.
This process is repeated until all the solutions are include in a front level.

###  4.4 Offspring generation

The offspring generation is related on how the individuals of the new
population are selected. In the case of NSGA-II, an elitism technique is used
konak2006multi , where the new individuals are only included in the new
generation if they improve the most dominant solutions.

In each algorithm iteration, a set of new individuals, with the same size than
the population, is generated (lines 11 to 17 in Algorithm 1). These new
individuals are joined with the individuals in the parent population (line
18), obtaining a set of the double size of the population. Only the half of
the individuals need to be selected for the offspring. The front levels and
crowding distances of the new set are calculated and used to order the
solution set (lines 19 to 22). The best half of the solutions is selected as
individuals of the offspring (line 23), and the rest are rejected.

###  4.5 Execution parametrization

Finally, the parametrization of the algorithm is required. The values of these
parameters also influence the performance of the algorithm. Previous testing
and/or bibliography research need to be done to establish them.

The parametrization is related to the starting and ending conditions of the
algorithm and how often some of the chromosome operators are applied. In our
case, the initial population is generated randomly (line 5 in Algorithm 1).
The number of individuals in the population is fixed to 200 solutions. The
ending condition is established by the number of generations. We fix this
value in 300 generations. Finally, the crossover probability is fixed to 1.0
because the design of NSGA-II already considers the possibility of keeping
fathers from one generation to the offspring. The mutation probability is
establish in 0.25 and the selection of the mutation is done uniformly. Table 2
summarizes the parameters for the execution of the genetic algorithm.

Table 2: Genetic algorithm execution parameters Parameter | Value
---|---
Population size | 200
Generations number | 300
Mutation probability | 0.25
Mutation selection probability | $1/3$
Crossover probability | $1.0$
Algorithm 1 Multi-objective optimization algorithm

1:procedure NSGA-II

2: $populationSize\leftarrow 200$

3: $generationNumber\leftarrow 300$

4: $mutationProb\leftarrow 0.25$

5: $P_{t}\leftarrow generateRandomPopulation(populationSize)$

6: $fitness\leftarrow calculateFitness(P_{t})$

7: $fronts\leftarrow calculateFronts(P_{t},fitness)$

8: $distances\leftarrow calculateCrowding(P_{t},fronts,fitness)$

9: for $i<generationNumber$ do

10: $P_{off}=\emptyset$

11: for $j<populationSize$ do

12: $father1\leftarrow binaryTournamentSelection(P_{t},fronts,distances)$

13: $father2\leftarrow binaryTournamentSelection(P_{t},fronts,distances)$

14: $child1,child2\leftarrow crossover(father1,father2)$

15: if $random()<mutationProb$ then $mutation(child1)$

16: if $random()<mutationProb$ then $mutation(child2)$

17: $P_{off}=P_{off}\cup\\{child1,child2\\}$

18: $P_{off}=P_{off}\cup P_{t}$

19: $fitness\leftarrow calculateFitness(P_{t})$

20: $fronts\leftarrow calculateFronts(P_{off},fitness)$

21: $distances\leftarrow calculateCrowding(P_{off},fronts,fitness)$

22: $P_{off}=orderElements(P_{off},fronts,distances)$

23: $P_{t}=P_{off}[1..populationSize]$ #the best half

24: $Solution=fronts[1]$ #the Pareto front

##  5 Experimental evaluation

In our experiments, we measured the values of the multi-objective functions
for different levels of workload and cluster capacity and the number of
generations needed to find suitable values. Our results were also compared
with another approach.

The first step for the execution of the experiments was to set up the values
for the parameters of our model. The access to real data from enterprise’s
applications is not alway straight forward. Thus, the parameters for the model
were estimated from scenarios as similar as possible to our application case.
Tables 3 and 4 show the parameters values for system and workload models and
for the microservices stack, respectively.

The selection of the parameters values for the microservices were based on the
analysis of Socks Shop weaveworks2016shocksshop . Socks Shop is a
microservices demo application that was developed to test the benefits of
deploying applications in a container platform. The demo also includes a load
test that we used to measure the computational resources consumption for each
microservice by benchmarking a live demo111http://cloud.weave.works/demo. The
values for the failure rates of the microservices were estimated using data
from cloud application traces di2015gloudsim  and the size of the source code
of Socks Shop222https://github.com/microservices-
demo333https://hub.docker.com/u/weaveworksdemos/. The number of microservices
requests for each application request was estimated from the customer behavior
model graph (CBMG) of an e-commerce benchmark menasce2002tpc .

We considered a heterogeneous cluster with four types of physical machines
that differed between them in their computational capacity: $cap=[100.0,$
$200.0,$ $400.0,$ $800.0]$. The four machines types were fixed with the same
value for the fail rate, $fail=0.025$ fu2010failure . The experiments
contained the same quantity of machines of each type. The physical machines
were distributed in two racks. The network distance for physical machines in
the same rack was fixed in $d=1.0$. Network distance for machines in different
racks was fixed in $d=4.0$.

In order to study our approach under different conditions, we executed the
optimization process for several experiment configurations. The workload of
the system was varied between the experiments by modifying the number of
applications running in the system ($|A|=[1,2]$) and by changing the number of
user requests that an application received ($ureq_{j}=[1.0,1.5,2.0]$). The
capacity of the cluster was modified by changing the number of available
physical machines ($|P|=[250,300,350,400]$). Thus, our experimentation
included 24 configurations.

Table 3: Experiments set up Element | Parameter | Value
---|---|---
Application | $|A|$ | 1 and 2
| $ureq_{j}$ | 1.0, 1.5 and 2.0
Physical machine | $|P|$ | 250, 300, 350 and 400
| $cap_{l}$ | 100.0, 200.0, 400.0 or 800.0
| $fail_{l}$ | 0.025
Network | $dist_{pm_{l},pm_{l^{\prime}}}$ | 1.0 or 4.0
Microservice |  | Table 4
Table 4: Microservices stack for Socks Shop application Name | Id. | Consumes | $msreq_{i}$ | $res_{i}$ | $thr_{i}$ | $fail_{i}$
---|---|---|---|---|---|---
worker | $ms_{1}$ | $\\{\\}$ | 3.2 | 0.1 | 1.0 | 0.04
shipping | $ms_{2}$ | $\\{ms_{13}\\}$ | 1.8 | 11.7 | 25.0 | 0.02
queue-master | $ms_{3}$ | $\\{s_{2},s_{4}\\}$ | 3.2 | 20.0 | 200.0 | 0.02
payment | $ms_{4}$ | $\\{\\}$ | 1.4 | 0.1 | 10.0 | 0.0002
orders | $ms_{5}$ | $\\{ms_{2},ms_{4},ms_{10},ms_{11},ms_{12},\\}$ | 2.3 | 27.1 | 80.0 | 0.02
login | $ms_{6}$ | $\\{\\}$ | 0.8 | 2.8 | 30.0 | 0.0001
front-end | $ms_{7}$ | $\\{ms_{5},ms_{6},ms_{9}\\}$ | 15.1 | 3.8 | 50.0 | 0.003
edge-router | $ms_{8}$ | $\\{ms_{7}\\}$ | 15.1 | 0.5 | 10.0 | 0.0001
catalogue | $ms_{9}$ | $\\{\\}$ | 12.0 | 0.2 | 3.0 | 0.0006
cart | $ms_{10}$ | $\\{ms_{12}\\}$ | 3.2 | 41.3 | 100.0 | 0.02
accounts | $ms_{11}$ | $\\{ms_{12}\\}$ | 0.1 | 45.1 | 100.0 | 0.003
weavedb | $ms_{12}$ | $\\{\\}$ | 3.2 | 26.3 | 80.0 | 0.04
rabbitmq | $ms_{13}$ | $\\{ms_{1},ms_{3}\\}$ | 3.2 | 4.0 | 40.0 | 0.0006
consul | $ms_{14}$ | $\\{\\}$ | 3.2 | 13.2 | 100.0 | 0.0003

We did not considered the similar researches explained in Section LABEL: to
compare our solution with, because their differences against our approach.
Therefore, the comparison of our results was done with Kubernetes container
policies.

Kubernetes uses two policies to manage the allocation of containers or pods,
as they are called in Kubernetes implementation: PodFitResources and
LeastRequestPriority vohra2017scheduling . The first one limits that the sum
of requested resources on a physical machine must not be greater than the
capacity of the machine. The second one determines that containers are
allocated in machines with least consumed resources. Both policies are focused
on spreading the containers across all the machines in the cluster. These
policies can not manage the scale level of the containers. In order to a
fairer comparison of the results, the scale level of the experiments for
Kubernetes was fixed at the same values that our approach obtained in the
experiment results. The other experiment configurations and conditions were
exactly the same than in the experiments for our approach.

##  6 Results

The results were analyzed from two perspectives: the adaptation of the
solutions to the objectives functions as new generations were obtained
(evolution of the objectives functions values) and the results of the
optimization —final values for the optimization objectives—. In the case of
the first perspective, a summary of three representative cases is presented.

It is important to remember that the solution obtained by a multi-objective
optimization evolutionary algorithm in general, and the NSGA-II in particular,
is not a single solution but a set of solutions that correspond to all the
individuals placed in the Pareto optimal front. In favor of simplifying the
analysis of the results, three data series are represented: (a) the minimum
value of an objective between all the solutions in the Pareto front, labeled
as min; (b) the mean of all the values of an objective for all the solutions
in the Pareto front, labeled as mean; and (c) the value of the objective for a
representative solution in the Pareto front, labeled as minSOV. The selection
of the representative solution was done by choosing the one with the smaller
single objective value (SOV). This SOV was calculated as the weighted average
of the normalized values of the four objectives (Equation 17). Our criteria
for selecting a solution was arbitrary and it was done only with the propose
of giving more details in the analysis of the results. The criteria to choose
between all the set of solutions with the same quality, the Pareto front, will
depend on the conditions and requirements of the specific environment where
the solution will be implemented. For instance, an environment in which
reliability of the system is more important than other objectives, solutions
in the Pareto front with smaller values for the System Failure objective
should be selected.

$\displaystyle
SOV=0.25\frac{value_{threshold}-min_{threshold}}{max_{threshold}-min_{threshold}}+0.25\frac{value_{cluster}-min_{cluster}}{max_{cluster}-min_{cluster}}$
(17)
$\displaystyle+0.25\frac{value_{failure}-min_{failure}}{max_{failure}-min_{failure}}+0.25\frac{value_{network}-min_{network}}{max_{network}-min_{network}}$
(18)

Figure 4: Evolution of the four objectives functions values along the
generations for the experiment with 1 application, with a workload level of
1.0 user requests using 300 physical machines

Figure 5: Evolution of the system state (nodes and services) along the
generations for the experiment with 1 application, with a workload level of
1.0 user requests using 300 physical machines

Figure 6: Evolution of the four objectives functions values along the
generations for the experiment with 2 application, with a workload level of
1.5 user requests using 350 physical machines

Figure 7: Evolution of the system state (nodes and services) along the
generations for the experiment with 2 application, with a workload level of
1.5 user requests using 350 physical machines

Figure 8: Evolution of the four objectives functions values along the
generations for the experiment with 1 application, with a workload level of
2.0 user requests using 400 physical machines

Figure 9: Evolution of the system state (nodes and services) along the
generations for the experiment with 1 application, with a workload level of
2.0 user requests using 400 physical machines

In the first case, the analysis of the evolution of the solutions, the results
for three experiment configurations are presented: Figures 4 and 5 for the
experiment with 300 physical machines, 1.0 user requests and 1 application;
Figures 6 and 7 for 350 physical machines, 1.5 user requests and 2
applications; and Figures 8 and 9 for 400 physical machines, 2.0 user requests
and 1 application.

Results are presented separately for the objectives functions —Figures 4, 6
and 8, including the series min, mean and minSOV— and for the system
conditions —Figures 5, 7 and 9, including the series min and mean—. System
conditions are represented with microservices scale level and number of
physical machines that allocate at least one microservice.

Finally, one of the previous configurations (350 physical machines, 1.5 user
requests and 2 application) is used to illustrate the evolution of all the
solutions in the Pareto set (Figures 10 and 11). We represent the scatter
plots of the values of the objectives functions for six generations. As each
of the points in the scatter plot is a four values tuple, one for each
objective function, two scatter plots are used for the representation: one for
the values of Total Network Distance and Cluster Balanced Used (Figure 10) and
another one for System Failure and Threshold Distance (Figure 11).

Figure 10: Evolution of the Pareto Optimal Front along the generations for
objectives Network Distance and Cluster Balanced Used for the experiment with
2 applications, 1.5 user requests and 350 physical machines

Figure 11: Evolution of the Pareto Optimal Front along the generations for
objectives System Failure and Threshold Distance for the experiment with 2
applications, 1.5 user requests and 350 physical machines

In the second case, the analysis of final optimization results, we present the
results for all the 24 experiment configurations. These results are plotted
separately for each optimization function and for the total number of used
physical machines: Figure 12 for the Total Network Distance; Figure 13 for
Threshold Distance; Figure 14 for the Cluster Balanced User; Figure 15 for the
System Failure; and Figure 16 for the number of nodes with allocated
containers. To represent a specific value for an experiment and an objective
function, one solution among the solutions included in the Pareto front needs
to be selected. We used again the criteria based on selecting the solution
with the smaller SOV fitness value as we have previously explained. The
solutions using our approach, labeled as NSGA-II, are compared with the
solutions obtained using the allocation policies of Kubernetes, labeled as
Kubernetes.

Figure 12: Comparative results for Network Distance optimization

Figure 13: Comparative results for Network Distance optimization

Figure 14: Comparative results for Cluster Balanced Used optimization

Figure 15: Comparative results for System Failure optimization

Figure 16: Comparative results for System Failure optimization

##  7 Discussion

The analysis of the results, as in the case of their presentation in the
previous section, is done from two points of views: the analysis of evolution
of the objectives and the analysis of final optimized results.

###  7.1 Evolution of the objectives

In general terms, it is observed that the values for the objectives functions
are established at the end of the execution of the optimization algorithm.
This stabilization is even achieved before generation number 100 (Figures 4, 6
and 8).

The min series shows a very fast stabilization for the System Failure,
Threshold Distance and Network Distance objectives. The min series for Cluster
Balanced is stabilized fairly later. In the case of the mean series, it is
observed that the values of the objectives are increased for some of them and
decreased for other ones before they meet a stabilization. The reason is that
the dominant solutions show equilibrated balances between all the objectives,
thus some of them need to be deteriorated to allow the other objectives to be
improved.

We want to highlight the behavior of the representative solution of the Pareto
front, the one with the smaller SOV value. In all the three cases, the value
of the solution minSOV is better than the values for the mean series for some
of the objectives and is worse for other objectives. For instance, in Figure 8
the minSOV solution is better for the System Failure, the Threshold Distance
and the Network Distance, but worse for the Cluster Balanced. This represents
again the fact that for improving one of the objectives, it is necessary to
degrade other objectives. It is important to remember that the values of the
min series for different objectives belong to different solutions and they
represent the lowest value obtained considering all the solutions of our
approach.

Finally, in relation with the number of physical machines with allocated
containers (Figures 5, 7 and 9), it is observed that the number of machines in
use is smaller than the available ones for all the configurations. This is
also a benefit, because other optimizations as energy consumption reduction
can be achieve by switching off the machines with no allocated containers.
Furthermore, the number of microservices represented in the figures indicates
that the scale level for each container is between 18 and 25 for the mean and
between 2 and 7 for the min, considering that there are 14 microservices for
each application. This high level of scalability is explained because it
benefits the reduction of the fail rate of the system, but it does not
necessarily damage the other objectives, except for the Threshold Distance. A
future improvement of the optimization process could be to include some type
of constraint to avoid solutions with too high levels of scalability.

Figures 10 and 11 show how the initial number of solutions in the Pareto front
is very small and how it increases as new generations are considered. We want
to highlight how just in 30 generations the solutions in the Pareto front are
shifted from random positions to the positions corresponding to the solution
search space. From this generation, the solution space is not shifted but the
number of solutions in the solution space is increased. This can be seen very
clearly in Figure 11, where due to the uniform distribution of the containers
in the initial population, all the solutions in the initial Pareto front are
concentrated in only one point of Threshold Distance and System Failure. In
just 30 generations, the Pareto front is covering all the solution space, and
in the following generations the solution space is populated with a bigger
diversity of solutions. By this, the diversity of solutions with the same
quality is increased and the covered solution space is bigger as new
generations are obtained.

We want to highlight that the solution search spaces of three of the
objectives are well covered. But in the case of the System Failure objective,
the solutions are not covering all the possible values (Figure 11). To
guarantee a better distribution of the solutions across the solution search
space in this objective, additional mutation operators should be considered.

In general terms, we can state that the use of a genetic approach based on
NSGA-II is suitable to find optimized solutions with a reasonable number of
generation (100) and a reasonable size of the population (200).

###  7.2 Final optimized results

In a second step, the final optimized results from all the experiment
configurations are analyzed. The results for our approach correspond to the
optimization values for the solution with a smaller SOV (Figures 12, 13, 14
and 15). In general terms, we can observe that for the objectives Network
Distance and Threshold Distance, the cases with the worst (higher)
optimization values correspond to the ones with a configuration of 2
applications: as the number of microservices is increased, these two
objectives values are worse. This behavior is also observed in the case of
increasing the workload of the system (number or user application requests).
On the contrary, the Cluster Balanced Use and the System Failure do no show a
clear trend as the workload and number of applications is varied. Finally, the
number of physical machines (Figure 16) seems to be influenced in the same way
that the two first objectives functions: the total number of physical nodes is
increased as the workload or the number of applications increase.

In relation with the results obtained with the Kubernetes approach, our
solution showed better optimizations for all the objectives. It is important
to highlight that these improvements were obtained in spite of using a smaller
number of physical machines. This also supposed a collateral benefit because
the machines with none allocated container could be switched off obtaining
greater reduction in the power consumption of the cluster than in the case of
Kubernetes.

##  8 Conclusion

In this article, we have addressed the problem of container resource
allocation and elasticity (self-scaling of container replicas). A four-
objectives optimization has been proposed based on a uniform distribution of
the workload along the physical machines (Cluster Balanced Use), a tight
distribution of the microservices workload along their container replicas
(Threshold Distance), a reduction of the network overheads (Network Distance)
and a better reliability (System Failure).

We have applied the NSGA-II, a genetic algorithm, to the problem of container
allocation in cloud. Several experiments configurations, varying the number of
physical machines and the workload of the system, have been tested in order to
study the benefits of this algorithm. The results have shown that our approach
is a suitable solution to address this problem and it found optimized
solutions with reasonable number of generations (100) and reasonable size of
the population (200).

The results have been compared with the obtained from the implementation of
the Kubernetes’ allocation policies. Our approach has shown a better behavior
by obtaining a better optimization of the four objectives. Moreover, our
solution uses a smaller number of physical machines.

As future work, we plan to study the results of our allocation algorithm in a
real Cloud Container Cluster. Additionally, other optimization objectives can
been included. Finally, container live-migration cost could be studied as the
criteria to select solutions from the Pareto front, instead of using a
criteria based on a single scalar value of the multi-objectives.

## References

  * (1) Resource management for infrastructure as a service (iaas) in cloud computing: A survey.  Journal of Network and Computer Applications 41, 424 – 440 (2014).  DOI http://dx.doi.org/10.1016/j.jnca.2013.10.004.  URL http://www.sciencedirect.com/science/article/pii/S1084804513002099
  * (2) Akhter, N., Othman, M.: Energy aware resource allocation of cloud data center: review and open issues.  Cluster Computing 19(3), 1163–1182 (2016)
  * (3) Balalaie, A., Heydarnoori, A., Jamshidi, P.: Microservices architecture enables devops: Migration to a cloud-native architecture.  IEEE Software 33(3), 42–52 (2016).  DOI 10.1109/MS.2016.64
  * (4) Beloglazov, A., Abawajy, J.H., Buyya, R.: Energy-aware resource allocation heuristics for efficient management of data centers for cloud computing.  Future Generation Comp. Syst. 28(5), 755–768 (2012).  DOI 10.1016/j.future.2011.04.017.  URL http://dx.doi.org/10.1016/j.future.2011.04.017
  * (5) Caballer, M., Blanquer, I., Moltó, G., de Alfonso, C.: Dynamic management of virtual infrastructures.  Journal of Grid Computing 13(1), 53–70 (2015).  DOI 10.1007/s10723-014-9296-5.  URL http://dx.doi.org/10.1007/s10723-014-9296-5
  * (6) Deb, K., Pratap, A., Agarwal, S., Meyarivan, T.: A fast and elitist multiobjective genetic algorithm: Nsga-ii.  Trans. Evol. Comp 6(2), 182–197 (2002).  DOI 10.1109/4235.996017.  URL http://dx.doi.org/10.1109/4235.996017
  * (7) Di, S., Cappello, F.: Gloudsim: Google trace based cloud simulator with virtual machines.  Softw. Pract. Exper. 45(11), 1571–1590 (2015).  DOI 10.1002/spe.2303.  URL http://dx.doi.org/10.1002/spe.2303
  * (8) Fazio, M., Celesti, A., Ranjan, R., Liu, C., Chen, L., Villari, M.: Open issues in scheduling microservices in the cloud.  IEEE Cloud Computing 3(5), 81–88 (2016).  DOI doi.ieeecomputersociety.org/10.1109/MCC.2016.112
  * (9) Fu, S.: Failure-aware resource management for high-availability computing clusters with distributed virtual machines.  Journal of Parallel and Distributed Computing 70(4), 384 – 393 (2010).  DOI http://dx.doi.org/10.1016/j.jpdc.2010.01.002.  URL http://www.sciencedirect.com/science/article/pii/S0743731510000031
  * (10) Guan, X., Wan, X., Choi, B.Y., Song, S., Zhu, J.: Application oriented dynamic resource allocation for data centers using docker containers.  IEEE Communications Letters PP(99), 1–1 (2016).  DOI 10.1109/LCOMM.2016.2644658
  * (11) Kang, D.K., Choi, G.B., Kim, S.H., Hwang, I.S., Youn, C.H.: Workload-aware resource management for energy efficient heterogeneous docker containers.  In: 2016 IEEE Region 10 Conference (TENCON), pp. 2428–2431 (2016).  DOI 10.1109/TENCON.2016.7848467
  * (12) Konak, A., Coit, D.W., Smith, A.E.: Multi-objective optimization using genetic algorithms: A tutorial.  Reliability Engineering and System Safety 91(9), 992 – 1007 (2006).  DOI http://dx.doi.org/10.1016/j.ress.2005.11.018.  URL http://www.sciencedirect.com/science/article/pii/S0951832005002012.  Special Issue - Genetic Algorithms and Reliability Special Issue
  * (13) von Lücken, C., Barán, B., Brizuela, C.: A survey on multi-objective evolutionary algorithms for many-objective problems.  Computational Optimization and Applications 58(3), 707–756 (2014).  DOI 10.1007/s10589-014-9644-1.  URL http://dx.doi.org/10.1007/s10589-014-9644-1
  * (14) Menasce, D.A.: Tpc-w: a benchmark for e-commerce.  IEEE Internet Computing 6(3), 83–87 (2002).  DOI 10.1109/MIC.2002.1003136
  * (15) Mitchell, M.: An Introduction to Genetic Algorithms.  MIT Press, Cambridge, MA, USA (1998)
  * (16) Pascual, J.A., Lorido-Botrán T.and Miguel-Alonso, J., Lozano, J.A.: Towards a greener cloud infrastructure management using optimized placement policies.  Journal of Grid Computing 13(3), 375–389 (2015).  DOI 10.1007/s10723-014-9312-9.  URL http://dx.doi.org/10.1007/s10723-014-9312-9
  * (17) Peinl, R., Holzschuher, F., Pfitzer, F.: Docker cluster management for the cloud - survey results and own solution.  Journal of Grid Computing 14(2), 265–282 (2016).  DOI 10.1007/s10723-016-9366-y.  URL http://dx.doi.org/10.1007/s10723-016-9366-y
  * (18) Singh, S., Chana, I.: Cloud resource provisioning: survey, status and future research directions.  Knowl. Inf. Syst. 49(3), 1005–1069 (2016).  DOI 10.1007/s10115-016-0922-3.  URL http://dx.doi.org/10.1007/s10115-016-0922-3
  * (19) Singh, S., Chana, I.: A survey on resource scheduling in cloud computing: Issues and challenges.  J. Grid Comput. 14(2), 217–264 (2016).  DOI 10.1007/s10723-015-9359-2.  URL http://dx.doi.org/10.1007/s10723-015-9359-2
  * (20) Vohra, D.: Scheduling pods on nodes.  In: Kubernetes Management Design Patterns, pp. 199–236. Springer (2017)
  * (21) Weaveworks, ContainerSolutions: Socks shop - a microservices demo application (2016).  URL https://microservices-demo.github.io/
  * (22) Wei, G., Vasilakos, A.V., Zheng, Y., Xiong, N.: A game-theoretic method of fair resource allocation for cloud computing services.  The Journal of Supercomputing 54(2), 252–269 (2010).  DOI 10.1007/s11227-009-0318-1.  URL http://dx.doi.org/10.1007/s11227-009-0318-1
  * (23) Xu, H., Liu, W., Shu, G., Li, J.: Location-aware data block allocation strategy for hdfs-based applications in the cloud.  2016 IEEE 9th International Conference on Cloud Computing 00, 252–259 (2016).  DOI doi.ieeecomputersociety.org/10.1109/CLOUD.2016.0042
  * (24) Xu, X., Yu, H., Pei, X.: A novel resource scheduling approach in container based clouds.  In: Computational Science and Engineering (CSE), 2014 IEEE 17th International Conference on, pp. 257–264. IEEE (2014)
  * (25) Zhan, Z.H., Liu, X.F., Gong, Y.J., Zhang, J., Chung, H.S.H., Li, Y.: Cloud computing resource scheduling and a survey of its evolutionary approaches.  ACM Comput. Surv. 47(4), 63:1–63:33 (2015).  DOI 10.1145/2788397.  URL http://doi.acm.org/10.1145/2788397
  * (26) Zhang, Q., Cheng, L., Boutaba, R.: Cloud computing: state-of-the-art and research challenges.  Journal of Internet Services and Applications 1(1), 7–18 (2010).  DOI 10.1007/s13174-010-0007-6.  URL http://dx.doi.org/10.1007/s13174-010-0007-6


<EMAIL_ADDRESS>
University of Würzburg

# A numerical solution approach for non-smooth optimal control problems based
on the Pontryagin maximum principle

Daniel Wachsmuth1  [

###### Abstract

We consider nonsmooth optimal control problems subject to a linear elliptic
partial differential equation with homogeneous Dirichlet boundary conditions.
It is well-known that local solutions satisfy the celebrated Pontryagin
maximum principle. In this note, we will investigate an optimization method
that is based on the maximum principle. We prove that the discrepancy in the
maximum principle vanishes along the resulting sequence of iterates. Numerical
experiments confirm the theoretical findings.

##  1 Introduction

In this note, we consider the following optimal control problem: Minimize

$J(y,u):=\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}+\int_{\Omega}g(u(x))\,\text{{}{}d}x$
(1)

over all $u\in L^{2}(\Omega)$ and $y\in H^{1}_{0}(\Omega)$ satisfying

$\displaystyle-\Delta y$ $\displaystyle=u$ $\displaystyle\text{ in }\Omega,$
$\displaystyle y$ $\displaystyle=0$ $\displaystyle\text{ on }\partial\Omega.$

Here, $\Omega\subseteq\mathbb{R}^{d}$ is a bounded domain, and
$g:\mathbb{R}\to\bar{\mathbb{R}}=\mathbb{R}\cup\\{+\infty\\}$ is assumed to be
proper and lower semicontinuous. In addition, we require

$\lim_{|v|\to\infty}\frac{g(v)}{|v|}=+\infty.$ (2)

Note, that we assume neither continuity nor convexity of $g$. Hence, it is
impossible to prove existence of solutions of (1). In fact, one can construct
problems without solution, see [18, Section 4.5]. In this note, we will work
with the example

$g(u):=\frac{\alpha}{2}u^{2}+I_{\mathbb{Z}}(u)=\begin{cases}\frac{\alpha}{2}u^{2}&\text{
if }u\in\mathbb{Z}\\\ +\infty&\text{ otherwise,}\end{cases}$ (3)

where $\alpha>0$. If $g$ is assumed to be convex and continuous, then
existence of solutions of (1) can be proven by the direct method of the
calculus of variations [17]. Let us remark that by the above assumptions $g$
is bounded from below.

If solutions exist, then the Pontryagin maximum principle [11] is a necessary
optimality condition. Its main feature is that no differentiability with
respect to the controls is needed, and so it is perfectly suited for the
problems considered here. In fact, due to the structure of the problem (linear
state equation, convexity of $J$ with respect to $y$), the maximum principle
is sufficient. We refer to [2, 5, 4, 12] for the Pontryagin maximum principle
applied to optimal control problems for partial differential equations. The
goal of this note is to construct an algorithm to solve the maximum principle.
We will comment on related work in Section 4.

##  2 Sensitivity analysis

In this section, we will perform a sensitivity analysis with respect to
perturbations of the control with characteristic functions. The setup is as
follows: Let $u,\tilde{u}\in L^{2}(\Omega)$ be feasible controls, i.e., the
integrals $\int_{\Omega}g(u)\,\text{{}{}d}x$ and
$\int_{\Omega}g(\tilde{u})\,\text{{}{}d}x$ exist. Let $B\subseteq\Omega$ be
measurable. We define

$u_{B}:=u+\chi_{B}(\tilde{u}-u).$

Let $y,y_{B}$ be the uniquely determined weak solutions of

$\displaystyle-\Delta y$ $\displaystyle=u$ $\displaystyle-\Delta y_{B}$
$\displaystyle=u_{B}$ $\displaystyle\text{ in }\Omega,$ $\displaystyle y$
$\displaystyle=0$ $\displaystyle y_{B}$ $\displaystyle=0$ $\displaystyle\text{
on }\partial\Omega.$

Let $p\in H^{1}_{0}(\Omega)$ be the weak solution of the adjoint equation

$\displaystyle-\Delta p$ $\displaystyle=y-y_{d}$ $\displaystyle\text{ in
}\Omega,$ $\displaystyle p$ $\displaystyle=0$ $\displaystyle\text{ on
}\partial\Omega.$

The goal is now to estimate $J(y_{B},u_{B})-J(y,u)$ in terms of
$u,\tilde{u},p$ and the Lebesgue measure $|B|$ of $B$. Here, we have the
following result.

###### Lemma 2.1

Under the assumptions above, we have

$J(y_{B},u_{B})-J(y,u)=\int_{B}(\tilde{u}-u)p+g(\tilde{u})-g(u)\,\text{{}{}d}x+\frac{1}{2}\|y_{B}-y\|_{L^{2}(\Omega)}^{2}.$

###### Proof 2.2

This follows directly from the definition of $p$ and $u_{B}$:

$\displaystyle J(y_{B},u_{B})-J(y,u)$
$\displaystyle=\frac{1}{2}\|y_{B}-y_{d}\|_{L^{2}(\Omega)}^{2}+\int_{\Omega}g(u_{B})\,\text{{}{}d}x-\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}-\int_{\Omega}g(u)\,\text{{}{}d}x$
$\displaystyle=\int_{\Omega}(y_{B}-y)(y-y_{d})+\frac{1}{2}(y_{B}-y)^{2}\,\text{{}{}d}x+\int_{B}g(\tilde{u})-g(u)\,\text{{}{}d}x$
$\displaystyle=\int_{B}(\tilde{u}-u)p+g(\tilde{u})-g(u)\,\text{{}{}d}x+\frac{1}{2}\|y_{B}-y\|_{L^{2}(\Omega)}^{2}.$

We will now prove that $\|y_{B}-y\|_{L^{2}(\Omega)}^{2}$ is of higher order
with respect to the Lebesgue measure $|B|$ of $B$.

###### Lemma 2.3

There are constants $c>0$ and $\nu>1/2$ independent of $u,\tilde{u},B$ such
that

$\|y_{B}-y\|_{L^{2}(\Omega)}\leq
c\,|B|^{\nu}\cdot\|\tilde{u}-u\|_{L^{\infty}(\Omega)},$

where $|B|$ denotes the Lebesgue measure of $B$. The constant $\nu$ can be
chosen as

$\nu=\begin{cases}1&\text{ if }d\leq 3,\\\ 1-\epsilon&\text{ if }d=4\text{ for
}\epsilon>0,\\\ \frac{1}{2}+\frac{2}{d}&\text{ if }d>4.\end{cases}$

###### Proof 2.4

We prove the claim by a well-known duality argument. Assume $d\leq 3$. Let
$w\in L^{2}(\Omega)$ be given. Let $z,q\in H^{1}_{0}(\Omega)$ be the weak
solutions of

$\displaystyle-\Delta z$ $\displaystyle=w$ $\displaystyle-\Delta q$
$\displaystyle=z$ $\displaystyle\text{ in }\Omega,$ $\displaystyle z$
$\displaystyle=0$ $\displaystyle q$ $\displaystyle=0$ $\displaystyle\text{ on
}\partial\Omega.$

Due to [15], there is $c>0$ independent of $w,z$ such that

$\|z\|_{L^{\infty}(\Omega)}\leq c\|w\|_{L^{2}(\Omega)}.$

Testing the weak formulations with $z$ and $q$ yields

$\|z\|_{L^{2}(\Omega)}^{2}=\int_{\Omega}wq\,\text{{}{}d}x\leq\|w\|_{L^{1}(\Omega)}\|q\|_{L^{\infty}(\Omega)}\leq
c\|w\|_{L^{1}(\Omega)}\|z\|_{L^{2}(\Omega)}.$

This proves $\|z\|_{L^{2}(\Omega)}\leq c\|w\|_{L^{1}(\Omega)}$. Applying this
estimate to $z:=y_{B}-y$ and $w:=u_{B}-u$ yields the claim with

$\|y_{B}-y\|_{L^{2}(\Omega)}\leq c\|u_{B}-u\|_{L^{1}(\Omega)}\leq
c\,|B|\cdot\|\tilde{u}-u\|_{L^{\infty}(\Omega)}.$

In case $d>3$ one can use the estimates from [3, Theorem 18].

Combining these results proves the following theorem.

###### Theorem 2.5

Let $u,\tilde{u}\in L^{\infty}(\Omega)$. Let $B\subseteq\Omega$ be measurable.
Let $\tilde{u},y_{B},y,p$ be defined as above. Then there are $\gamma>0$ and
$c>0$ independent of $u,\tilde{u},B$ such that

$J(y_{B},u_{B})-J(y,u)\leq\int_{B}(\tilde{u}-u)p+g(\tilde{u})-g(u)\,\text{{}{}d}x+c\,|B|^{1+\gamma}\|\tilde{u}-u\|_{L^{\infty}(\Omega)}^{2}.$

##  3 Pontryagin maximum principle

With the help of Theorem 2.5 we can prove the Pontryagin maximum principle.

###### Theorem 3.1

Let $\bar{u}\in L^{\infty}(\Omega)$ be locally optimal with respect to
$L^{1}(\Omega)$ topology for the control problem (1). Let $\bar{y},\bar{p}\in
H^{1}_{0}(\Omega)$ be the optimal state and adjoint solving

$\displaystyle-\Delta\bar{y}$ $\displaystyle=\bar{u}$
$\displaystyle-\Delta\bar{p}$ $\displaystyle=\bar{y}-y_{d}$
$\displaystyle\text{ in }\Omega,$ $\displaystyle\bar{y}$ $\displaystyle=0$
$\displaystyle\bar{p}$ $\displaystyle=0$ $\displaystyle\text{ on
}\partial\Omega.$

Let $v\in\mathbb{\mathbb{R}}$ be such that $g(v)<+\infty$. Then

$\bar{u}(x)\bar{p}(x)+g(\bar{u}(x))\leq v\bar{p}(x)+g(v)\text{ for almost all
}x\in\Omega.$ (4)

###### Proof 3.2

Let $v\in\mathbb{\mathbb{R}}$ be such that $g(v)<+\infty$. Applying Theorem
2.5 with $u:=\bar{u}$, $\tilde{u}:=v$ yields

$0\leq
J(y_{B},u_{B})-J(\bar{y},\bar{u})=\int_{B}(v-\bar{u})\bar{p}+g(\tilde{u})-g(\bar{u})\,\text{{}{}d}x+o(|B|).$

By standard arguments based on the Lebesgue differentiation theorem, see,
e.g., [10, Theorem 2.1], the claim follows.

The maximum principle is a sufficient condition for the problem considered
here.

###### Corollary 3.3

Let $\bar{u}\in L^{2}(\Omega)$ satisfy the conclusion (4) of Theorem 3.1. Then
$\bar{u}$ is global optimal for (1).

###### Proof 3.4

Let $\tilde{u}\in L^{2}(\Omega)$ be an admissible control with associated
state $\tilde{y}$. Then Lemma 2.1 with $B=\Omega$ yields

$J(\tilde{y},\tilde{u})-J(\bar{y},\bar{u})=\int_{\Omega}(\tilde{u}-\bar{u})\bar{p}+g(\tilde{u})-g(\bar{u})\,\text{{}{}d}x+\frac{1}{2}\|\tilde{y}-\bar{y}\|_{L^{2}(\Omega)}^{2}.$

Since $\bar{u}$ satisfies (4), the first expression is non-negative, which
implies
$J(\tilde{y},\tilde{u})-J(\bar{y},\bar{u})\geq\frac{1}{2}\|\tilde{y}-\bar{y}\|_{L^{2}(\Omega)}^{2}\geq
0$.

##  4 Construction of an algorithm

We will now apply Theorem 2.5 with $u:=u_{k}$ and $\tilde{u}:=\tilde{u}_{k}$,
where $u_{k}$ is the current iterate of the algorithm to be devised. Let
$y_{k}$ and $p_{k}$ be the associated state and adjoint. The control
$\tilde{u}_{k}$ has to be computed in each iteration. Let $B_{k}$ be
measurable. Then we have

$J(y_{B_{k}},u_{B_{k}})-J(y_{k},u_{k})=\int_{B_{k}}(\tilde{u}_{k}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x+o(|B_{k}|).$
(5)

The idea is now to choose $\tilde{u}_{k}$ and $B_{k}$ such that
$J(y_{B_{k}},u_{B_{k}})-J(y_{k},u_{k})$ is negative and to define the new
iterate by

$u_{k+1}=u_{k}+\chi_{B_{k}}(\tilde{u}_{k}-u_{k}).$

In view of the maximum principle, Theorem 3.1, it is natural to choose
$\tilde{u}_{k}$ as a function satisfying

$\tilde{u}_{k}(x)\in\operatorname*{arg\,min}_{v\in\mathbb{R}}vp_{k}+g(v).$ (6)

In addition, $B_{k}$ will be chosen to get sufficient descent.

Let us comment on related work. The classic algorithm of [8] chooses
$B_{k}:=\Omega$, resulting in a fixed-point scheme to solve the maximum
principle. The min-h method of [7] uses the update
$u_{k+1}:=u_{k}+t(\tilde{u}_{k}-u_{k})$ with $t\in(0,1]$, and is thus only
suited for convex functions $g$. In the monograph [14], a method similar to
ours is presented to solve optimal control problems with ODEs. Let us also
also mention the review papers [6, 16]. In [9] binary control problems are
solved with a similar approach: there a trust-region globalization is
proposed, whereas we use an Armijo line-search to globalize.

As motivated above, we will compute $\tilde{u}_{k}$ as a result of the
pointwise minimization

$\tilde{u}_{k}(x)\in\operatorname*{arg\,min}_{v\in\mathbb{R}}vp_{k}+g(v).$

Due to (2) this problem is solvable for all $x$. A measurable selection of
this argmin-map exists [1]. For the example of $g$ proposed in (3), we get

$\tilde{u}_{k}(x)\in\operatorname{round}\left(-\frac{1}{\alpha}p_{k}(x)\right).$

It remains to describe how $B_{k}$ is chosen. Here, we are faced with two
competing goals: In order to make the first term in (5) as small as possible,
$B_{k}$ has to be chosen as large as possible. However, to control the
remainder term in (5), $|B_{k}|$ has to be chosen sufficiently small.

We propose the following line-search. Given $t\in(0,1]$, choose $B_{t}$ such
that

$\displaystyle\int_{B_{t}}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x$
$\displaystyle\leq
t\int_{\Omega}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x,$
(7) $\displaystyle|B_{t}|$ $\displaystyle\leq t\cdot|\Omega|.$

Due to the celebrated Lyapunov convexity theorem, see, e.g., [13, Theorem
5.5], a measurable set $B_{t}$ satisfying (7) exists. Given $t$ and $B_{t}$,
we set $u_{t}:=u_{k}+\chi_{B_{t}}(\tilde{u}_{k}-u_{k})$. Let $y_{t}$ be the
associated state.

The parameter $t_{k}$ is determined by the following procedure: Let $t_{k}$ be
the largest number in $\\{\beta^{l}:\ l\in\mathbb{N}\cup\\{0\\}\\}$, where
$\beta\in(0,1)$, that satisfies the descent condition

$J(y_{t},u_{t})-J(y_{k},u_{k})\leq\sigma\int_{B_{t}}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x$
(8)

where $\sigma\in(0,1)$, and $B_{t}$ is a measurable set satisfying (7). This
condition is inspired by the well-known Armijo line-search in nonlinear
optimization. If $u_{k}$ does not satisfy the maximum principle, there is an
admissible step-size $t_{k}$, and the resulting algorithm produces a new
iterate with smaller value of the objective.

###### Lemma 4.1

Suppose that

$\int_{\Omega}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x<0.$

There is $t_{0}>0$ such that for all $t\in(0,t_{0})$ condition (8) is
satisfied.

###### Proof 4.2

Due to Theorem 2.5, we have

$J(y_{t},u_{t})-J(y_{k},u_{k})-\sigma\int_{B_{t}}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x\\\
\leq(1-\sigma)\int_{B_{t}}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x+o(t)\\\
\leq
t(1-\sigma)\int_{\Omega}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x+o(t),$

which proves the claim.

The resulting algorithm is sketched in Algorithm 1.

Choose $\beta\in(0,1)$, $\sigma\in(0,1)$, $u_{0}$ with
$\int_{\Omega}g(u_{0})\,\text{{}{}d}x<\infty$, $\delta_{\mathrm{tol}}\geq 0$.
Set $k:=0$.

loop$\triangleright$ Gradient descent

Compute state $y_{k}$ and adjoint $p_{k}$ associated to $u_{k}$.

Compute $\tilde{u}_{k}$ as in (6).

if
$\left|\int_{\Omega}(\tilde{u}_{k}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x\right|\leq\delta_{\mathrm{tol}}$
then $\triangleright$ Termination criterion

return $u_{k}$

end if

$t:=1$.

loop$\triangleright$ Armijo line-search

Compute $B_{k,t}$ satisfying (7).

Compute $J(y_{t},u_{t})$.

if (8) is satisfied then

break

end if

$t:=\beta\cdot t$.

end loop

$t_{k}:=t$. $\triangleright$ Update

$u_{k+1}:=u_{k}+\chi_{B_{k,t_{k}}}(\tilde{u}_{k}-u_{k})$.

$k:=k+1$.

end loop

Algorithm 1 Maximum-principle based descent algorithm

Let us now turn to the convergence analysis of Algorithm 1. Here, we follow
the related analysis in [19]. Let us define

$\rho_{k}:=\int_{\Omega}(\tilde{u}_{k}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x.$

Due to the choice of $\tilde{u}_{k}$ in (6), it follows $\rho_{k}\leq 0$. If
$\rho_{k}=0$ then $u_{k}$ satisfies the maximum principle Theorem 3.1, and the
corresponding control $u_{k}$ is optimal by Corollary 3.3.

###### Lemma 4.3

Let $(u_{k})$ be an infinite sequence generated by Algorithm 1. Then

$\sum_{k=0}^{\infty}t_{k}\|\rho_{k}\|_{L^{1}(\Omega)}<+\infty.$

###### Proof 4.4

Using conditions (8) and (7) shows

$J(y_{k+1},u_{k+1})-J(y_{k},u_{k})\leq\sigma\int_{B_{t_{k}}}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x\leq
t_{k}\int_{\Omega}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x=-t_{k}\|\rho_{k}\|_{L^{1}(\Omega)}.$

Due to (2), $g$ has a global minimum and is bounded from below, so that $J$ is
bounded from below by some $M\in\mathbb{R}$. Summing this inequality over
$k\in\mathbb{N}$ and using $J\geq M$ proves
$\sum_{k=1}^{\infty}t_{k}\|\rho_{k}\|_{L^{1}(\Omega)}\leq
J(y_{0},u_{0})-M<\infty$.

For simplicity, we assume for the subsequent convergence analysis that

$\operatorname{dom}g:=\\{v:\ g(v)<\infty\\}$ (9)

is compact. Then the set of iterates $(u_{k})$ and $(\tilde{u}_{k})$ is
uniformly bounded in $L^{\infty}(\Omega)$.

###### Corollary 4.5

Assume (9). Let $M>0$ such that $\operatorname{dom}g\subseteq[-M,+M]$. Then
$\|u_{k}\|_{L^{\infty}(\Omega)}\leq M$ and
$\|\tilde{u}_{k}\|_{L^{\infty}(\Omega)}\leq M$ for all $k$.

###### Theorem 4.6

Assume (9). Either the Algorithm 1 stops after finitely many steps with

$\left|\int_{\Omega}(\tilde{u}_{k}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x\right|\leq\delta_{\mathrm{tol}}$

(so that $u_{k}$ satisfies the maximum principle if
$\delta_{\mathrm{tol}}=0$), or

$\int_{\Omega}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x\to
0,$

i.e., the residual in the maximum principle tends to zero, and $(u_{k})$ is a
minimizing sequence.

###### Proof 4.7

We follow the proof of the related result [19, Theorem 6.7]. Let us assume the
algorithm generates an infinite sequence of iterates. Let $k$ be such that
$t_{k}<1$. Due to the line-search procedure of Algorithm 1, it follows that
$t:=\beta^{-1}t_{k}\leq 1$ violates the descent condition (8), that is

$0<J(y_{t},u_{t})-J(y_{k},u_{k})-\sigma\int_{B_{t}}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x.$

As in the proof of Lemma 4.1, we get from Theorem 2.5

$0<t(1-\sigma)\int_{\Omega}(\tilde{u}-u_{k})p_{k}+g(\tilde{u})-g(u_{k})\,\text{{}{}d}x+c\,|t|^{1+\gamma}\|\tilde{u}-u\|_{L^{\infty}(\Omega)}.$

Together with Corollary 4.5, we get

$0<-t(1-\sigma)\|\rho_{k}\|_{L^{1}(\Omega)}+c|t|^{1+\gamma},$

where $c$ is independent of $k$. This implies

$\|\rho_{k}\|_{L^{1}(\Omega)}\leq ct_{k}^{\gamma}$

for all $k$ such that $t_{k}<1$. With Lemma 4.3, we get

$+\infty>\sum_{k=0}^{\infty}t_{k}\|\rho_{k}\|_{L^{1}(\Omega)}=\left(\sum_{k:\,t_{k}=1}\|\rho_{k}\|_{L^{1}(\Omega)}\right)+\left(\sum_{k:\,t_{k}<1}t_{k}\|\rho_{k}\|_{L^{1}(\Omega)}\right)\geq\left(\sum_{k:\,t_{k}=1}\|\rho_{k}\|_{L^{1}(\Omega)}\right)+c\left(\sum_{k:\,t_{k}<1}\|\rho_{k}\|_{L^{1}(\Omega)}^{1+\frac{1}{\gamma}}\right),$

which results in $\lim_{k\to\infty}\|\rho_{k}\|_{L^{1}(\Omega)}=0$. Hence, the
algorithm stops after finitely many iterations if $\delta_{\mathrm{tol}}>0$.

##  5 Numerical results

Let us now report about numerical experiments with Algorithm 1. Here, we
consider the optimal control problem

$J(y,u):=\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}+\frac{\alpha}{2}\|u\|_{L^{2}(\Omega)}^{2}+I_{\mathbb{Z}\cap[-b,b]}(u)$

over all $u\in L^{2}(\Omega)$ and $y\in H^{1}_{0}(\Omega)$ satisfying

$\displaystyle-\Delta y$ $\displaystyle=u$ $\displaystyle\text{ in }\Omega,$
$\displaystyle y$ $\displaystyle=0$ $\displaystyle\text{ on }\partial\Omega.$

This fits into the setting of the paper with the choice

$g(v):=\frac{\alpha}{2}v^{2}+I_{\mathbb{Z}\cap[-b,b]}$

Here, we chose $\Omega=(0,1)^{2}$,

$y_{d}(x_{1},x_{2})=10x_{1}\sin(5x_{1})\cos(7x_{2}),\quad\alpha=0.01,\quad\beta=0.01,\quad
b=10.$

We discretized the problem with piecewise linear finite elements on a regular
mesh for state and adjoint variables, while the control was discretized with
piecewise constant finite elements. We report the results for a sequence of
different meshes, where the finest mesh has mesh-size $h=1.41\cdot 10^{-3}$
resulting in $\approx 2\cdot 10^{6}$ degrees of freedom for the control
variables, which results in a mixed-integer optimization problem with $\approx
2\cdot 10^{6}$ integer variables. In the implementation of Algorithm 1 a
greedy strategy was used to determine $B_{t}$. The loop in Algorithm 1 was
terminated if in the inner loop $t|\Omega|$ was smaller than any of the
elements in the grid.

Now let us report about some of the results. The optimal control can be seen
in the left plot of Figure 1. In the right plot, we report about the iteration
history of the residual $\|\rho_{k}\|_{L^{1}(\Omega)}$. Surprisingly, the
iterations seem to be mesh independent. In addition, for this particular
problem a very small number of iterations was needed to optimize over $2\cdot
10^{6}$ discrete control variables.

Figure 1: Optimal control (left), iteration history (right)

This is underlined by the results in Table 1. It shows for different
discretizations the final value of the objective $J$ and the final value of
the residual $\|\rho\|_{L^{1}(\Omega)}$. As can be seen from the last column
of this table, very few outer iterations are needed. In conclusion, this new
algorithm seems to be capable of solving quite challenging mixed-inter
programs.

$h$ | $J$ | $\|\rho\|_{L^{1}(\Omega)}$ | It
---|---|---|---
$4.42\text{\cdot}{10}^{-02}$ | $4.706$ | $3.20\text{\cdot}{10}^{-06}$ | $4$
$2.21\text{\cdot}{10}^{-02}$ | $5.048$ | $2.02\text{\cdot}{10}^{-08}$ | $6$
$1.13\text{\cdot}{10}^{-02}$ | $5.210$ | $6.00\text{\cdot}{10}^{-11}$ | $8$
$5.66\text{\cdot}{10}^{-03}$ | $5.293$ | $8.91\text{\cdot}{10}^{-11}$ | $8$
$2.83\text{\cdot}{10}^{-03}$ | $5.334$ | $6.46\text{\cdot}{10}^{-12}$ | $9$
$1.41\text{\cdot}{10}^{-03}$ | $5.354$ | $4.11\text{\cdot}{10}^{-13}$ | $10$
Table 1: Iteration history

## Acknowledgements

This research was supported by the German Research Foundation (DFG) under
grant number WA 3626/3-2 within the priority program “Non-smooth and
Complementarity-based Distributed Parameter Systems: Simulation and
Hierarchical Optimization” (SPP 1962). The author thanks Anna Lentz for
comments on an earlier version of this manuscript.

## References

  * [1] J.-P. Aubin and H. Frankowska.  Set-Valued Analysis, volume 2 of Systems & Control: Foundations & Applications.  Birkhäuser Boston, Inc., Boston, MA, 1990.  doi:10.1007/978-0-8176-4848-0.
  * [2] F. Bonnans and E. Casas.  An extension of Pontryagin’s principle for state-constrained optimal control of semilinear elliptic equations and variational inequalities.  SIAM J. Control Optim., 33(1):274–298, 1995.  doi:10.1137/S0363012992237777.
  * [3] H. Brézis and W. A. Strauss.  Semi-linear second-order elliptic equations in $L^{1}$.  J. Math. Soc. Japan, 25:565–590, 1973.  doi:10.2969/jmsj/02540565.
  * [4] E. Casas, J.-P. Raymond, and H. Zidani.  Pontryagin’s principle for local solutions of control problems with mixed control-state constraints.  SIAM J. Control Optim., 39(4):1182–1203, 2000.  doi:10.1137/S0363012998345627.
  * [5] E. Casas.  Pontryagin’s principle for optimal control problems governed by semilinear elliptic equations.  In Control and estimation of distributed parameter systems: nonlinear phenomena (Vorau, 1993), volume 118 of Internat. Ser. Numer. Math., pages 97–114. Birkhäuser, Basel, 1994.  doi:10.1007/978-3-0348-8530-0_6.
  * [6] F. L. Chernous’ko and A. A. Lyubushin.  Method of successive approximations for solution of optimal control problems.  Optimal Control Appl. Methods, 3(2):101–114, 1982.  doi:10.1002/oca.4660030201.
  * [7] R. G. Gottlieb.  Rapid convergence to optimum solutions using a min-h strategy.  AIAA Journal, 5(2):322–329, 1967.  doi:10.2514/3.3960.
  * [8] I. A. Krylov and F. L. Černous’ko.  The method of successive approximations for solving optimal control problems.  Ž. Vyčisl. Mat i Mat. Fiz., 2:1132–1139, 1962.  (In Russian).
  * [9] P. Manns, M. Hahn, C. Kirches, S. Leyffer, and S. Sager.  On convergence of binary trust-region steepest descent.  J. Nonsmooth Anal. Optim., 4, 2023.  doi:10.46298/jnsao-2023-10164.
  * [10] C. Natemeyer and D. Wachsmuth.  A proximal gradient method for control problems with non-smooth and non-convex control cost.  Comput. Optim. Appl., 80(2):639–677, 2021.  doi:10.1007/s10589-021-00308-0.
  * [11] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko.  The mathematical theory of optimal processes.  Interscience Publishers John Wiley & Sons, Inc., New York-London, 1962\.  Translated from the Russian by K. N. Trirogoff; edited by L. W. Neustadt.
  * [12] J. P. Raymond and H. Zidani.  Pontryagin’s principle for state-constrained control problems governed by parabolic equations with unbounded controls.  SIAM J. Control Optim., 36(6):1853–1879, 1998.  doi:10.1137/S0363012996302470.
  * [13] W. Rudin.  Functional analysis.  McGraw-Hill Series in Higher Mathematics. McGraw-Hill Book Co., New York-Düsseldorf-Johannesburg, 1973.
  * [14] V. A. Srochko.  Iterative methods for solving optimal control problems.  Fismatlit, Moscow, 2000.  (In Russian).
  * [15] G. Stampacchia.  Le problème de Dirichlet pour les équations elliptiques du second ordre à coefficients discontinus.  Ann. Inst. Fourier (Grenoble), 15(fasc. 1):189–258, 1965.  URL: http://www.numdam.org/item?id=AIF_1965__15_1_189_0.
  * [16] A. S. Strekalovsky.  Modern methods for solving nonconvex optimal control problems.  The Bulletin of Irkutsk State University. Series Mathematics, 8:141–163, 2014.  (In Russian).
  * [17] F. Tröltzsch.  Optimal control of partial differential equations, volume 112 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2010.  Theory, methods and applications, Translated from the 2005 German original by Jürgen Sprekels.  doi:10.1090/gsm/112.
  * [18] D. Wachsmuth.  Iterative hard-thresholding applied to optimal control problems with $L^{0}(\Omega)$ control cost.  SIAM J. Control Optim., 57(2):854–879, 2019.  doi:10.1137/18M1194602.
  * [19] D. Wachsmuth.  A topological derivative-based algorithm to solve optimal control problems with $L^{0}(\Omega)$ control cost.  J. Nonsmooth Anal. Optim., 5, 2024.  doi:10.46298/jnsao-2024-12366.