# A Preliminary Study for Identification of Additive Manufactured Objects with
Transmitted Images

Kenta Yamamoto<EMAIL_ADDRESS>University of
Tsukuba  Ryota Kawamura  University of Tsukuba  Kazuki Takazawa  Pixie Dust
Technologies, Inc.  Hiroyuki Osone  University of Tsukuba  Yoichi Ochiai
University of Tsukuba  Pixie Dust Technologies, Inc.

###### Abstract

Additive manufacturing has the potential to become a standard method for
manufacturing products, and product information is indispensable for the item
distribution system. While most products are given barcodes to the exterior
surfaces, research on embedding barcodes inside products is underway. This is
because additive manufacturing makes it possible to carry out manufacturing
and information adding at the same time, and embedding information inside does
not impair the exterior appearance of the product. However, products that have
not been embedded information can not be identified, and embedded information
can not be rewritten later. In this study, we have developed a product
identification system that does not require embedding barcodes inside. This
system uses a transmission image of the product which contains information of
each product such as different inner support structures and manufacturing
errors. We have shown through experiments that if datasets of transmission
images are available, objects can be identified with an accuracy of over 90%.
This result suggests that our approach can be useful for identifying objects
without embedded information.

##  1 Introduction

The additive manufacturing (AM) process has greatly improved the ease of model
customization. When managing products with a digital system, the addition of
information is essential, such as a barcode. While management information has
only to be attached to the exterior, additive manufacturing process makes it
possible to embed information inside the product. Detection is then possible
by various methods, for example, embedded information can be read using
terahertz (Thz) sensing in InfraStructs[1], and using visible light projection
in AirCode[2].

Two problems exist with the embedded information approach. The first is that
an information adding step must be introduced into the manufacturing process.
This also constitutes an advantage as the integration of manufacturing and
adding information enables overall reduction of processing time and cost,
nevertheless it is not possible to identify objects that were manufactured
without embedded information. The second problem arises when the embedded
information is duplicated such that product batches cannot be uniquely
identified. This can occur in the current item distribution system. In
previous research, defects in the barcode surface produced by the
manufacturing process have been used for unique identification of products[3].

In this study, a system of product identification was developed that does not
use embedded information. All products have inherent structure combined with
characteristics created by the manufacturing process, particularly
manufacturing errors. Transmission images of product items are utilized to
extract inherent product characteristics. This idea was inspired by research
that uses the finger vein pattern to authenticate individuals. Finger vein
pattern authentication is biometric, hence information is not added in advance
of the authentication process; and, in order to visualize the internal blood
vessels, an image of the finger transmitted by infrared light is used. Also in
our system, information is not added, and we use images transmitted by
infrared light.

To verify the identification system, two experiments were conducted: feature
matching and deep learning; using two types of object group with identical
exterior appearance. In one group, the internal structure of the object was
varied. In the other group, all structures, including the internal structure,
were unchanged; i.e., the objects differed only as a result of manufacturing
errors.

The results suggest that feature matching is possible when there is a change
in shape, and that deep learning can be used to identify objects that vary due
to naturally occurring manufacturing errors. The purpose of this investigation
is to examine the possibility of utilizing an identification system based on
transmission images, to discriminate products in the absence of embedded
information.

Figure 1: System Overview. The system is a three step process: fabrication, capturing, and identification. Products are fabricated without embedding information internally. Image capture requires transmitted light in the near infrared region. Identification is performed using feature matching or deep learning with datasets. Table 1: Comparison with existing methods. | Information Place | Pre Embedding | Detection Angle | Learning Dataset | Wavelength | Reflection or Transmission
---|---|---|---|---|---|---
InfraStructs[1] | Inside | Yes | One way | No | THz | Reflection
AirCode[2] | Inside | Yes | One way | No | VIS | Reflection
LayerCode[4] | Exterior | Yes | All angle | No | VIS or NIR | Reflection
Our System | Everywhere | No | All angle | Yes | NIR | Transmission

##  2 Related Works

###  2.1 Properties of Additive Manufacturing

Additive Manufacturing encompasses a variety of methods[5]. Fused Deposition
Modeling (FDM), in particular, has made a great contribution to the uptake of
the desktop 3D printer. The widespread adoption of FDM has been accompanied by
examination of the characteristics of FDM manufactured products. Many previous
studies have considered how the characteristics of objects are affected by the
manufacturing method itself; e.g., product characteristics[6, 7, 8, 9],
mechanical properties[10, 11], and material properties[12, 13, 14, 15].
Several investigations have also been conducted to determine how objects
change with different parameters; i.e., the relationship between each
parameter and mechanical properties[16, 17, 18, 19], and the correspondence
between the internal structure and mechanical properties[20, 21, 22, 23, 24].
By referring to this body of literature concerning the FDM method, the
fabrication characteristics of FDM could be established without the necessity
to carry out exploratory research on the process.

###  2.2 Biometrics using Transmitted Images

It has become very common to perform personal authentication using biometric
characteristics such as fingerprint authentication and face recognition. Most
biometric authentication methods use an image produced by reflection of light
from a target. However, some biometrics acquire an image after light has
passed through the target. Finger vein pattern recognition is a typical
example. In 2004, Miura et al.[25] published a paper on personal
authentication using a transmission image of a finger vein pattern obtained by
illumination with infrared LEDs. The main challenge for this method lies in
ensuring that the system is robust. It must be able to authenticate properly
even if finger placement changes slightly, and it must not authenticate as
another person. These problems led to research on how to automatically
determine the region of interest[26, 27, 28] and improve imaging quality[29,
30]. Therefore, research literature concerning transparent biometric
authentication was referred to in order to assist with the construction of an
identification system using transmission images of AM products.

###  2.3 Embedding Information with Additive Manufacturing

Additive Manufacturing is a promising rapid prototyping method since it
enables easy customization of shapes, thus prompting research on embedding
barcodes during manufacturing. InfraStructs[1] is a famous early study that
embodies this concept. During the manufacturing process, a structure was built
inside the object that was detected and decoded as information, using Thz
sensing. This information could not be captured with a visible light camera.
Several subsequent studies successfully embedded and read information in the
object interior using alternative methods: i.e., internal information read in
the far infrared region[31], information embedded using highly reflective
material[32], and information embedded in the form of air pockets and read
using visible light projection[2]. Recently, Maia et al.[4] embedded
information in the entire surface of an object in order to simplify reading of
the information.

Thus, existing research is based on building a system that embeds information
in an invisible form during manufacturing and reads it out. In contrast, a
method based on transmitted images does not require design of embedding
information and the reading angle is free. Table 1 summarizes these
differences. The utility of the transmission imaging system for identifying a
product without preprocessing has been verified.

Table 2: Fabrication characteristics. Main Topic |  Sub Topic |  Sub-Sub Topic
---|---|---
|  Printing Width |
|  Layer Thickness |
Shape |  |  Infill Pattern
Parameters |  Inner Support Structure |  Infill Density
|  |  Infill Position
|  |  Material
|  Material Properties |  Color
|  Printing Speed |
Condition |  Temperature of Heat-Bed |
Parameters |  Temperature of Printing Head |
|  Retraction Speed |
|  Environment |

##  3 Materials and Methods

###  3.1 System Overview

This section details the system overview. The system consists of three
processes: fabrication, capturing images, and identification (Fig.1). During
fabrication, information is not intentionally embedded in products. Instead,
inherent product characteristics are used for identification; i.e., visibly
indistinguishable products may differ in their internal structure, and
manufactured items hold intrinsic information such as dimensional values,
which differ due to manufacturing errors. In the capturing process,
transmission images are acquired. By using near infrared light, internal shape
is captured in the transmitted image to enable identification of products on
the basis of changes to internal structure. Inclusion of internal information
also improves the accuracy of identifying changes resulting from manufacturing
errors. Two methods were considered for the identification process: feature
matching and deep learning. Feature matching does not require pre-learning for
identification, allowing a simple system to be utilized. This contrasts with
deep learning, which does require pre-learning, however discrimination
accuracy is higher.

###  3.2 Fabrication Characteristics

Table 2 lists the characteristics that can appear in the transmission image of
the product. All characteristics are referred to in the existing white
paper[6] and the software that runs the 3D printer. The characteristics were
divided into two types, those related to shape and those related to condition.
Of the shape-related properties, we paid particular attention to the inner
support structure (see section 3.2.1). Many previous studies intentionally
modified the interior space of a product without affecting the exterior
appearance. In contrast, this study did not include intentional redesign of
the interior space, but instead focused on random internal changes that occur
naturally in manufacturing. Condition-related parameters were fixed in this
study as it was difficult to distinguish between the effects of changing these
parameters and differences arising from manufacturing errors. The assumption
was made that if a manufacturing error could be identified, a change in
condition parameters could also be distinguished. Therefore, only
manufacturing errors were considered.

Figure 2: Shape parameter variations used in the identification tests. (a) to
(c) Different infill patterns. (d) to (f) Different infill densities. (g) and
(h) Show the two in infill positions. (g) Initial infill position. (h)
Alternative infill position.

Inner Support Structure: This parameter consists of three sub-parameters:
infill pattern, infill density, and infill position (Fig.2), which are all
visible on the transmission image. Infill pattern is an important parameter
that is related to the shape of the internal structure and affects the
mechanical properties of the product. Modifying the infill pattern changes the
durability of the product under stress; i.e., the resistance to stress and the
strength of the product are modified. Infill density refers to the density of
a given infill pattern. Optimization of infill density is always desirable in
terms of manufacturing cost and cycle time, as the amount of material
increases with density. Infill position refers to the placement of the inner
support structure for the given infill pattern. As placement is often done
automatically by the slicer software, in general the system user is less aware
of infill position.

###  3.3 Transparent Imaging System

The transmissive imaging system operates by capturing the characteristics and
internal structure of products, both primitive and complex, in images using
light wavelengths that are transmitted through the target object.

####  3.3.1 Optical Setup

Figure 3: Optical System. (a) A picture of the imaging system used in this
study. (b) A plan view of (a), showing the path of light. (c) Alternative
imaging system.

An optical system using near infrared (NIR) light was constructed to obtain an
optimal transmission image of the object. The target object, NIR camera, NIR
light source, and an optical element were arranged as shown in Fig.3 (b). In
order to eliminate the influence of speckle, the halogen lamp was used as the
near infrared light source, which was equipped with an infrared transmission
filter. A Dihedral Corner Reflector Array (DCRA)[33] was used to irradiate the
back of the object. As shown in Fig.3 (b), DCRA has the function of imaging
the light source at a plane-symmetrical position. Similar functions can be
realized with mirror and lens as shown in Fig.3 (c). In this experiment,
because it was difficult to incorporate large mirror and lens, the amount of
illuminating light was maximized by using DCRA.

####  3.3.2 Transmitted Images

In Section 3.2, it was explained that inner support structures differ in type,
and that for a given support structure, the infill patterns, densities, and
positions of structural elements are variable. In this section, the appearance
of the transmission image is examined when the inner support structure is
either changed by design, or left unchanged such that variation arises from
manufacturing errors.

Infill Pattern: Cube and bunny shapes were manufactured using three types of
infill pattern: Diamond Fill, Hexagonal, and Linear. Transmission image
results are shown in Fig.4. A different image was obtained for each infill
pattern, for both cube and bunny.

Figure 4: Variation of additive manufacturing infill pattern. (top) The
transmission image when the infill pattern of the cube is changed. (bottom)
The transmission image when the infill pattern of the bunny is changed. Figure
5: Variation of additive manufacturing infill density. (top) The transmission
image when the infill density of the cube is changed. (bottom) The
transmission image when the infill density of the bunny is changed.

Infill Density: Cube and bunny shapes were produced with Diamond Fill pattern
densities of 5, 10, and 15% for the cube, and 10, 20, and 30% for the bunny.
As the density increased, the amount of transmitted light decreased, hence the
contrast in the transmissive images also decreased. These results are shown in
Fig.5.

Infill Position: The same inner structure was created at multiple positions,
with a Diamond Fill inner pattern density of 10% for the cube, and 20% for the
bunny. Transmission image results are shown Fig.6. The cube images change as
the position of the inner support structure changes. However, the transmission
images of the bunny were almost unchanged due to the complexity of the model.

Figure 6: Variation of additive manufacturing infill position (top) The
transmission image when the infill position of the cube is changed. (bottom)
The transmission image when the infill position of the bunny is changed.
Figure 7: Random variation in items made with the same manufacturing
parameters. (top) The comparison of transmission images of cubes with the same
infill conditions. (bottom) The comparison of transmission images of additive
manufactured bunnies with the same infill conditions.

Same Inner Parameters (only Manufacturing Errors): We compared the
transmission images when all the parameters are equal. In order to identify
objects manufactured with the same parameters, it is necessary to utilize
slight changes such as manufacturing errors. We observed whether such a small
difference could be acquired within a transmission image. As shown in Fig.7,
we prepared the cubes and bunnies with all same parameters listed in Table 2.
Since very similar images are acquired, the region of interest is enlarged. In
the enlarged image, it can be confirmed that trivial differences occur in both
the cube and the bunny.

##  4 Results

Figure 8: The results of feature matching when: (a) infill pattern is changed;
(b) infill density is changed; (c) infill position is changed; and (d)
comparing three bunnies with the same infill condition.

In Section 3.3, the relationships between the characteristics of objects and
their transmission images were established. When parameters related to shape
or structure differ between objects, clear differences are likely to appear
between their transmission images. However, only minor differences are visibly
discernible between transmission images of objects with the same parameters.
To evaluate the level of difference that can be resolved between various
images, two algorithm-based experiments were conducted: feature matching, and
classification by deep learning.

###  4.1 Feature Matching

To verify whether or not differences between the NIR transmission images were
detectable by image processing, feature extraction and feature matching were
initially conducted. In this experiment, the SURF algorithm[34] was used to
perform feature extraction.

The results of matching feature points are shown in Fig.8. Three key features
were extracted for use in feature point matching. The matching features in a
reference image and a target image were obtained by the k-nearest neighbor
method. The D. Lowe ratio test[35] was used to select the top two matched
features, the remaining feature is drawn in Fig.8. Feature matching was then
performed between pairs of images that varied by infill pattern, or infill
density, or infill position, or only by random differences due to
manufacturing errors. Of the paired images showing matching results, the image
on the left is the reference image, and the image on the right is the target
image. The right hand images were changed and the results of feature matching
were recorded as the number of matching points after the ratio test. These
results are listed in Table 3. They quantify the change in matching with the
change in a parameter; i.e., the degree to which the algorithm recognises the
difference between images. A single reference image is used for all variations
of a parameter hence the locations and maximum number of feature matching
points are retained in the parameter series, thus allowing quantification and
direct comparison of the results.

The results in Fig.8 and Table 3 demonstrate that it is possible to
discriminate between objects with obvious shape or structural differences by
feature extraction and feature point matching using the SURF algorithm. This
is not the case when objects are manufactured with the same structural
parameters and differ only by manufacturing errors. Quantitative differences
in matching results from one target image to the next were not significant,
therefore feature matching is unlikely to be able to distinguish between
objects at a higher degree of similarity.

Table 3: Match rate of feature points. Since the same image is used to compare the same parameters, the match rate is 100%. Infill Pattern |  Diamond |  Linear |  Hexagonal
---|---|---|---
Diamond |  255 / 255 |  82 / 255 |  62 / 266
Linear |  79 / 261 |  261 / 261 |  86 / 261
Hexagonal |  67 / 305 |  81 / 305 |  305 / 305
Infill Density |  10% |  20% |  30%
10% |  310 / 310 |  17 / 310 |  3 / 310
20% |  11 / 255 |  255 / 255 |  3 / 255
30% |  2 / 174 |  5 / 174 |  174 / 174
Infill Position |  Position A |  Position B |  Position C
Position A |  258 / 258 |  58 / 258 |  54 / 258
Position B |  42 / 215 |  215 / 215 |  35 / 215
Position C |  49 / 265 |  39 / 265 |  265 / 265
Same Position |  Object 1 |  Object 2 |  Object 3
Object 1 |  235 / 235 |  128 / 235 |  107 / 235
Object 2 |  122 / 227 |  227 / 227 |  78 / 227
Object 3 |  104 / 204 |  88 / 204 |  204 / 204

###  4.2 Deep Learning

Figure 9: Flow of dataset creation for deep learning. (top) How to move an object when creating a dataset. (bottom) An example of a transmission image that was learned. Table 4: Dataset Type. | 3D Printer | Position | Pattern | Density | Amount
---|---|---|---|---|---
a | Makerbot | Random | Diamond Fill | 20% | 10
b | Makerbot | Fixed | Diamond Fill | 20% | 10
c | Makerbot | Fixed | Diamond Fill | 10% | 10
d | Bellulo | Fixed | Honeycomb | 10% | 10
Figure 10: Deep learning second test. (top) Graph of learning progress when
there is no noise in the dataset. (bottom) Graph of learning progress when
noise is added to the dataset.

Deep learning has contributed greatly to the improvement of classification
accuracy in recent years. Following the partial success of feature matching as
discussed in the previous section, the deep learning technique was adopted to
distinguish between items that are made using the same parameters, but differ
due to manufacturing errors. Four groups of products were manufactured and
their transmission images acquired. Datasets were produced from the
transmission images (Table 4). Dataset (a) consists of products with different
infill positions. Datasets (b), (c), and (d) consist of products that differ
only as a result of manufacturing errors.

In order to apply the deep learning method, training datasets must first be
prepared. The imaging system described in Section 3.3 (Fig. 3) was used to
obtain images for the datasets. In a commercial or industrial setting, the
target may rotate slightly, or the set position can shift. To ensure that the
identification system would be sufficiently robust to cope with this type of
image variation, images with the target rotated and translated were included
in the training dataset. As shown in Fig.9, 25 points were set at 0.5 mm
intervals. By moving to each point and rotating 360 degrees at that location,
images of objects transmitted at all angles were acquired. Transmission images
were also obtained continuously, using the camera in video mode while the
target was moving and rotating. Datasets were created by extracting an image
from each frame of the video data.

Deep learning identification tests were performed using the datasets described
above. Resnet[36] was employed as the learning model and was trained in
advance to improve learning efficiency.

Datasets for the first learning test consisted only of the data from standard
transmission images; i.e., no added rotation or noise. The progress of the
accuracy rate is shown in Fig.10(left). The accuracy of discrimination between
similar images was high from the outset, due to the use of a trained neural
network. As the classification rate reached 100% in some cases, a more
difficult classification task was set for the next test.

For the second learning test, the dataset included pre-processed images
(augmented data), e.g., Gaussian noise and rotation images (5, 10, and 15
degrees). The progress of the accuracy rate is shown in Fig.10(right).
Although the level of accuracy was lower than achieved in the first learning
test, it was sufficient to distinguish between similar images. Hence,
identification by deep learning is effective, even with noise added to the
dataset.

##  5 Discussions

The proposed system was able to identify products with high accuracy in the
identification trials, however, there are still many points to be discussed.
The limitations of the study are described below.

###  5.1 LIMITATIONS AND IMPROVEMENTS

Transmitted Light and Contrast: The experimental investigation described in
this paper was conducted on the premise that the transmitted light and
contrast resolution, as available from the optical setup, were sufficient to
acquire the characteristics of the objects being imaged. In future, it will be
necessary to determine the minimum required transmitted light intensity, and
the threshold of the contrast required to resolve all details in transmitted
images.

Product Shape: As a preliminary investigation, this study was limited, on the
basis of time and cost, to using a cube as a simple model and a bunny as a
complex model. Further verification of the deep learning method will require
resources for a large number of products to be manufactured with specified
conditions, and their transmission images to be acquired as datasets. Products
that do not have a flat bottom will also require a fixture to enable NIR
photography.

Unverified Parameters: Changes to the “Inner Support Structure” sub-
parameters, and manufacturing errors were the only product modifications
considered in this paper. While differences arising from both were verified as
recognizable by the identification system, the effects of changing the
remaining unverified parameters must also be evaluated. “Printing Width” and
“Layer Thickness” are additional shape-related parameters (Table 2). Neither
was tested in this investigation as both can affect the overall appearance of
the product. However, in some product designs it is possible to change these
parameters without affecting the appearance of the product. The effects of
varying the parameters related to condition are also still to be tested.

Product Aging: In a real world distribution system, the product surface may be
scratched or otherwise damaged during transportation, and product color can
change over time. Although internal structure is less exposed to external
impacts, the inner shape and the inner materials can undergo aging effects.
The level of response of the detection system to such changes must be
evaluated.

Optical System: The optical system is shown in Fig.3. Although DCRA was used
to increase the amount of backlight in the small scale optical setup, the
optical system should not require such special optical elements. In this
study, the DCRA often appeared in transmission images (Fig.8). This problem
needs to be addressed since transmission images are used as training data for
the neural network, i.e. transmitted data must consist only of the object
image, without any form of extraneous information. The optical system requires
optimization in accordance with these issues.

Object Quantity: As shown in Table 4, 10 objects were prepared for each
dataset. This is because it took a considerable amount of time from the
production of objects to shooting videos for datasets, hence it was difficult
to prepare a large number of objects. However, the actual item distribution
system often handles more than 1000 items. In the future, it is necessary to
study how to deal with such a large amount of data.

Deep Learning Downside: Deep learning was adopted to realize high-precision
identification. This method always creates the necessity for preparation of
datasets and pre-learning. Dataset creation is also inefficient within the
current system. The process requires simplification.

Figure 11: This is a future application that, like a point of sales (POS)
system, identifies an object in real time and returns the result when the
object is placed.

##  6 Conclusion

A product identification system using transmission images was developed.
Product parameters relevant to the system were defined, and identification
experiments were conducted based on a subset of several selected parameters.
An optical system for acquiring transmission images was constructed, and the
acquired transmission images were observed. Two methods, feature matching and
deep learning, were verified as identification methods, and their
functionality and limitations were analyzed.

Future experiments will be conducted under more complex conditions, using
parameters not yet verified and objects of various shapes. Improvements will
be made to the system, enabling it to operate effectively under many
conditions. It will then be a useful system for identifying additively
manufactured products to which the information embedding method cannot be
applied (Fig.11).

## References

  * [1] Karl D. D. Willis and Andrew D. Wilson.  Infrastructs: fabricating information inside physical objects for imaging in the terahertz region.  ACM Trans. Graph., 32:138:1–138:10, 2013.
  * [2] Dingzeyu Li, Avinash S. Nair, Shree K. Nayar, and Changxi Zheng.  Aircode: Unobtrusive physical tags for digital fabrication.  ArXiv, abs/1707.05754, 2017.
  * [3] Rina Ueno and Jin Mitsugi.  Barcode fingerprinting: Unique identification of commercial products with their jan/ean/ucc barcode.  pages 416–420, 02 2018.
  * [4] Henrique Teles Maia, Dingzeyu Li, Yuan Yang, and Changxi Zheng.  Layercode: optical barcodes for 3d printed shapes.  ACM Trans. Graph., 38:112:1–112:14, 2019.
  * [5] Kaufui V Wong and Aldo Hernandez.  A review of additive manufacturing.  ISRN Mechanical Engineering, 2012, 2012.
  * [6] T. Grimm.  Fused deposition modeling: A technology evaluation.  Time-Compression Technologies, pages 1–6, 01 2003.
  * [7] Nur Saaidah Abu Bakar, Mohd Rizal Alkahari, and Hambali Boejang.  Analysis on fused deposition modelling performance.  Journal of Zhejiang University-SCIENCE A, 11:972–977, 2010.
  * [8] P. Roshan Kumar, Inderpreet Singh Ahuja, and Tejbir Singh.  Application of fusion deposition modelling for rapid investment casting – a review.  2012\.
  * [9] David A. Roberson, David Espalin, and Ryan B. Wicker.  3d printer selection: A decision-making evaluation and ranking model.  2013\.
  * [10] Antonio Lanzotti, Marzio Grasso, Gabriele Staiano, and M. Martorelli.  The impact of process parameters on mechanical properties of parts fabricated in pla with an open-source 3-d printer.  2015\.
  * [11] Jonathan Torres, Matthew Cole, Allen Owji, Zachary DeMastry, and Ali P. Gordon.  An approach for mechanical property optimization of fused deposition modeling with polylactic acid via design of experiments.  2016\.
  * [12] Ludmila Novakova-Marcincinova and Ivan Kuric.  Basic and advanced materials for fused deposition modeling rapid prototyping technology.  2012\.
  * [13] Jonathon Schofield and Michael Dawson.  Evaluation of dimensional accuracy and material properties of the makerbot 3d desktop printer.  Rapid Prototyping Journal, 21:618 – 627, 11 2015.
  * [14] Sunpreet Singh, Seeram Ramakrishna, and Rupinder Singh.  Material issues in additive manufacturing: A review.  Journal of Manufacturing Processes, 25:185–200, 01 2017.
  * [15] Jason Cantrell, Sean Rohde, David Damiani, Rishi Gurnani, Luke DiSandro, Josh Anton, Andie Young, Alex Jerez, Douglas Steinbach, Calvin Kroese, and Peter Ifju.  Experimental characterization of the mechanical properties of 3d-printed abs and polycarbonate parts.  Rapid Prototyping Journal, 23, 06 2017.
  * [16] Angel R. Torrado and David A. Roberson.  Failure analysis and anisotropy evaluation of 3d-printed tensile test specimens of different geometries and print raster patterns.  Journal of Failure Analysis and Prevention, 16:154–164, 2016.
  * [17] Behzad Rankouhi, Sina Javadpour, Fereidoon Delfanian, and Todd Letcher.  Failure analysis and mechanical characterization of 3d printed abs with respect to layer thickness and orientation.  Journal of Failure Analysis and Prevention, 16:467–481, 2016.
  * [18] J.M. Chacón, Miguel Caminero, Eustaquio García-Plaza, and Pedro López.  Additive manufacturing of pla structures using fused deposition modelling: Effect of process parameters on mechanical properties and their optimal selection.  Materials & design, 124:143–157, 03 2017.
  * [19] Diana Popescu, Aurelian Zapciu, Catalin Gheorghe Amza, Florin Baciu, and Rodica Marinescu.  Fdm process parameters influence over the mechanical properties of polymer specimens: A review.  Polymer Testing, 69, 05 2018.
  * [20] Liseli Jeanette Baich, Guha Manogharan, and Hazel Marie.  Study of infill print design on production cost-time of 3d printed abs parts.  2016\.
  * [21] Fernandez-VicenteMiguel, CalleWilson, FerrandizSantiago, and ConejeroAndres.  Effect of infill parameters on tensile mechanical behavior in desktop 3d printing.  2016\.
  * [22] C KennyLAlvarez, C RodrigoFLagos, and Miguel Aizpun.  Investigating the influence of infill percentage on the mechanical properties of fused deposition modelled abs parts.  2016\.
  * [23] Can Mert Al and Ulas Yaman.  Improving the strength of additively manufactured objects via modified interior structure.  2017\.
  * [24] Sridher Rangisetty and Larry Peel.  The effect of infill patterns and annealing on mechanical properties of additively manufactured thermoplastic composites.  page V001T08A017, 09 2017.
  * [25] Naoto Miura, Akio Nagasaka, and Takafumi Miyatake.  Feature extraction of finger-vein patterns based on repeated line tracking and its application to personal identification.  Machine Vision and Applications, 15(4):194–203, Oct 2004.
  * [26] Naoto Miura, Akio Nagasaka, and Takafumi Miyatake.  Extraction of finger-vein patterns using maximum curvature points in image profiles.  IEICE TRANSACTIONS on Information and Systems, 90(8):1185–1194, 2007.
  * [27] Wonseok Song, Taejeong Kim, Hee Chan Kim, Joon Hwan Choi, Hyoun-Joong Kong, and Seung-Rae Lee.  A finger-vein verification system using mean curvature.  Pattern Recognition Letters, 32:1541–1547, 2011.
  * [28] Jinfeng Yang and Yihua Shi.  Finger-vein roi localization and vein ridge enhancement.  Pattern Recognition Letters, 33:1569–1579, 2012.
  * [29] Kenichiro Tanaka, Yasuhiro Mukaigawa, Yasuyuki Matsushita, and Yasushi Yagi.  Descattering of transmissive observation using parallel high-frequency illumination.  IEEE International Conference on Computational Photography (ICCP), pages 1–8, 2013.
  * [30] Wenming Yang, Xiaola Huang, Fei Zhou, and Qingmin Liao.  Comparative competitive coding for personal identification by using finger vein and finger dorsal texture fusion.  Inf. Sci., 268:20–32, 2014.
  * [31] Ayumi Okada, Piyarat Silapasuphakornwong, Masahiro Suzuki, Hideyuki Torii, Youichi Takashima, and Kazutake Uehira.  Non-destructively reading out information embedded inside real objects by using far-infrared light.  In SPIE Optical Engineering + Applications, 2015.
  * [32] Masahiro Suzuki, Tomohisa Matumoto, Youichi Takashima, Hideyuki Torii, and Kazutake Uehira.  Information hiding inside 3-d printed objects by forming high reflectance projections.  In ICVIP, 2017.
  * [33] S Maekawa, K Nitta, and O Matoba.  Transmissive mirror device with micro dihedral corner reflector array-the mirror forming a real image.  ITE Journal, 30:49–52, 2006.
  * [34] Herbert Bay, Tinne Tuytelaars, and Luc Van Gool.  Surf: Speeded up robust features.  In Aleš Leonardis, Horst Bischof, and Axel Pinz, editors, Computer Vision – ECCV 2006, pages 404–417, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg.
  * [35] David G. Lowe.  Distinctive image features from scale-invariant keypoints.  International Journal of Computer Vision, 60(2):91–110, Nov 2004\.
  * [36] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2016.

Figure 12: Changes in transmission image due to layer characteristics. (a) a
photograph showing the experimental setup for imaging. (b) the actually
acquired transmission image. (c) an image of (b) displayed in pseudo color.
(d)-(g) the change of the transmission image when ”Layer Thickness” is
changed. (h), (i) the change of the transmission image when ”Printing Speed”
is changed.

##  Appendix A A Layer Feature

Figure 13: (left) a picture of an imaging system to verify robustness. (right
top) the imaging results when the XY stage is moved in X and Y directions.
(right middle) the imaging results when rotating the rotation stage. (right
bottom) the result when the camera is moved in X direction and imaged toward
the center of the object.

The transmission image of the product is a photograph of the light intensity
emitted from a certain point after passing through various paths. The smallest
unit of a product that constitutes such a phenomenon can be defined as a
”layer”. Regardless of the shape, infill pattern, infill density, and infill
position, the transmitted light has only passed through many layers. In other
words, the stacked feature of each layer becomes a transmission image. We have
already defined the characteristics of each layer in Section 3.2. Therefore,
here we verified how the transmission image changes according to the change in
the characteristics by changing the parameters of each layer.

The elements that make up the layer are broadly divided into two types: those
that affect the shape of the product and those that affect the state of the
product. The elements that affect the shape are ”Layer Thickness” and
”Printing Width”. The elements that affect the state are ”Printing Speed”,
”Temperature of Heat-Bed”, ”Temperature of Printing Head”, ”Retraction
Distance”, and ”Environment”. The reason why ”Inner Support Structure” is not
included in these elements is that ”Inner Support Strucure” is established by
the positional relationship of multiple layers.

###  A.1 Experimental Setup

In order to confirm the effect on the transmission image of one layer, we
installed a near infrared light source, a thin plate equivalent to one layer,
and a near infrared camera as shown in Fig. 12(a). To make it easier to see
how the light source changed before and after passing through the layer, the
shape of the light source entering the thin plate was made circular and the
light source was collimated. However, since coherent light is used as a light
source, the acquired image is accompanied by speckles. Therefore, in order to
give priority to the measurement of the rough change of the light source
shape, the influence of speckle was ignored and the saturation was allowed.
Also, in order to observe how the light diffuses after passing through one
layer, it is necessary to take a transmission image after passing through
multiple layers. In addition, we prepared a foundation with grooves formed at
intervals of 5 mm, and obtained a transmission image by fitting multiple
sheets into the foundation while changing the distance between the sheets.

###  A.2 Effect of Parameters

Of the two types of elements that make up the layer, changes in shape are
thought to have a large effect on the transmitted image. Therefore, we first
conducted an experiment to observe changes in the transmission image when the
parameters about the shape were changed. While ”Printing Width” has the same
value when using the same manufacturing equipment, ”Layer Thickness” can be
easily changed from the settings provided in the user interface. The
transmission images were compared. Fig.12(b) shows what changes appear in the
transmission image when the layer thickness is changed. In this experiment,
four layer thicknesses of 0.1 mm, 0.2 mm, 0.3 mm, and 0.4 mm were prepared. We
confirmed that the transmission image changed according to the thickness and
positional relationship of the layer.

Next, an experiment was conducted to observe the change in the transmitted
image accompanying the change in the parameters affecting the state. ”Printing
Speed” was selected as an element related to the state. The reason for
verifying only one of the many state parameters is that the change in the
state is considered to have a lower influence on the transmission image than
the change in the shape. Fig.12(h) and (i) show the imaging results. Compared
to when the shape parameters were changed (Fig.12(d)-(g)), the effect on the
transmitted image was small.

As described above, we have compared the transmission images of elements
related to shape and elements related to state. Although there was a
considerable change in both, the shape change is considered to be easy to use
as a clue for identification.

##  Appendix B Robust Transmissive Images

When using a transmission image for object identification, it is desirable
that the transmission image is not affected even if the orientation of the
backlight or the position and orientation of the object changes. In order to
achieve such robustness of the transmitted image, it is recommended that the
transmitted light sufficiently passes through multipath. In general, reflected
light, internal diffusion, and transmitted light are generated when light
enters an object. If sufficient internal diffusion occurs, the transmission
image obtained should be robust to the orientation and position of the target.

In order to verify the robustness of this system, we combined an XY scanning
stage and a rotating stage to change the orientation and position of the
target (Fig.13(left)). The XY scanning stage is Thorlabs MLS203-1, and the
rotary stage is Thorlabs PRMTZ8. In addition, for the object to be imaged, an
object with a primitive shape was prepared: a cube with a size of 50x50x50
[mm], an infill density of 10%, and an infill pattern of diamond fill.
Fig.13(right) show the imaging results when the XY scanning stage is moved,
the target is rotated, or the camera position is changed. Looking at the
region of interest in each imaging result, it was found that the orientation
and position of the target had little effect on the transmitted image.


soft open fences

# Piecewise linear value function approximations in nonlinear dynamic
scheduling problems with VTOLs

Viktoriya Nikitina
<EMAIL_ADDRESS>
&Sergejs Rogovs
<EMAIL_ADDRESS>
&Matthias Gerdts
<EMAIL_ADDRESS>

(Institute of Applied Mathematics and Scientific Computing
Department of Aerospace Engineering
University of the Bundeswehr Munich, Germany
October 29, 2022)

###### Abstract

Modern vertical take-off and landing vehicles (VTOLs) could significantly
affect the future of mobility. The broad range of their application fields
encompasses urban air mobility, transportation and logistics as well as
reconnaissance and observation missions in a military context. This article
presents an efficient numerical framework for dynamic scheduling of multiple
VTOLs. It combines a scheduling problem with optimal trajectory planning. The
whole setting is formulated as a mixed-integer bilevel optimization problem.
At the upper level, VTOLs are scheduled, and their starting times are
computed. The solution of the lower level problem involves the computation of
a value function and yields optimal trajectories for every aerial vehicle. In
order to solve the bilevel problem, it is recast into a single-level one. The
resulting mixed-integer nonlinear program (MINLP) is then piecewise linearized
and solved numerically by a linear solver based on the Branch-and-Bound
algorithm. Numerical results prove the feasibility of the approach proposed in
this work.

_Keywords_ Bilevel optimization problem $\cdot$ dynamic scheduling of VTOLs
$\cdot$ mixed-integer problem $\cdot$ piecewise linearization $\cdot$ shortest
path $\cdot$ value function.

##  1 Introduction

In the last few years, novel concepts of vertical take-off and landing
vehicles (VTOLs) have been developed and successfully tested in practice.
Modern approaches are considered to be promising and could considerably affect
the future of mobility. The recent technological advances in this field have
attracted much attention across different countries and industries. As a
result, the global VTOL market has lately seen a significant boom. In Europe,
for example, a number of requests have been recently received by the European
Union Aviation Safety Agency (EASA) for certifying new types of VTOLs [1].

Thanks to their high maneuverability, comparatively small size and universal
deployability (no need for runways), VTOLs are believed to play an important
role in the future of several industries. One of them is clearly urban air
mobility, where flying taxis could transport passengers from A to B in a much
faster way than ground vehicles due to the avoidance of congested streets.
Electric VTOLs (eVTOLs) could in addition yield environmental benefits and
promote sustainability in the aviation sector. One of the countries blazing a
trail in the field of air taxis is Singapore, whose government plans to start
the operation of eVTOLs by 2024 at one of the airports [2]. Another major area
of VTOL application relates to emergency situations, where aerial vehicles
could be used for quick medical evacuation of people injured in accidents or
for different search-and-rescue scenarios (e.g. after earthquakes, avalanches
etc.). Moreover, VTOLs could help to revolutionize the logistics sector. The
use of autonomous delivery drones could facilitate the whole shipping process
making it more efficient and cost-effective. Multiple key players in the field
of electronic commerce have already started to plan the integration of such
drones into their transport chain. In a military context, VTOLs could be, for
instance, involved in missions aimed at reconnaissance and observation of some
targets.

A rapidly increasing number of VTOLs requires their efficient coordination and
safe air traffic management, which might be particularly challenging when
dealing with limited space. This paper proposes an approach for a small fleet
of autonomous VTOLs (typically, up to eight) designed for their simultaneous
scheduling and optimal trajectory planning. The task of each VTOL is hereby to
reach a static target region without colliding with dynamic obstacles.
Provided that the starting position of each aerial vehicle is fixed and given,
the overall task is said to be completed as soon as the last VTOL has reached
the target set. The whole setting is formulated as a mixed-integer bilevel
optimization problem. At the upper level, the scheduling of multiple
autonomous VTOLs can be viewed as a job-shop problem (JSP) with the aim of
finding an optimal sequence and the corresponding starting times. At the lower
level, the task is to determine the shortest path from every aerial vehicle to
the target region. The solution procedure involves computation of a value
function and yields an optimal trajectory for every VTOL provided that the
starting times are given. Clearly, the optimization problems at both levels
are strongly coupled with each other and have to be therefore considered as a
whole.

To date, different efficient approaches have been developed to solve the
optimization problems at the upper and lower level if they are considered
separately. The authors of [3] and [4] deal, for instance, with optimal
trajectory planning for aircraft using the Hamilton-Jacobi-Bellman (HJB)
equation. [5] and [6] in turn present exact methods for coordination of
multiple landing aircraft. However, the class of coupled scheduling and
optimal control problems has not been extensively researched in the literature
to the best of the authors’ knowledge. In a context of autonomous vehicles,
[7] and [8] propose interesting methods for optimization problems with a
bilevel structure. Approaches presented in [9] and [10] are designed to
schedule multiple robots and provide an optimal path for each of them. In
[11], a bilevel mixed-integer optimization problem is investigated with the
aim of scheduling aircraft and planning their optimal trajectory in a busy
terminal area.

The remainder of this article is structured as follows. Section 2 briefly
introduces the value function of an optimal control problem, which solves the
Hamilton-Jacobi-Bellman (HJB) equation. This theoretical background originates
from dynamic programming and is required in succeeding sections. The
mathematical model used for dynamic scheduling of VTOLs is formulated in
Section 3. Its solution procedure is then carefully described in Section 4.
Herein, one of the steps includes an elegant technique used in this article
for the piecewise linearization of the final MINLP. The solution of the
resulting linearized problem can be obtained numerically using the robust and
efficient GUROBI solver. Section 5 is mainly devoted to the presentation of
the results achieved by the approach presented in this work. Finally, some
concluding remarks and possible directions for further research are given in
Section 6.

##  2 Optimal Control and Value Function

The main goal of this section is to give a short overview of the theory used
in Sections 3 and 4. The concepts briefly outlined below deal with the
Hamilton-Jacobi-Bellman (HJB) equation and the value function approach. While
Section 2.1 discusses the continuous time problem, Section 2.2 presents
formulas related to the discrete case. Note that the derivation of formulas is
omitted here since it goes beyond the scope of this article. More detailed
information can be, for example, found in [12, 13].

###  2.1 Continuous Time

For $\lambda\geq 0$ being a given constant discount rate,
$U\subset\mathbb{R}^{n_{u}}$ a control space (which is typically closed and
bounded), $x\in
W^{1,\infty}\left(\interval[openright]{0}{\infty},\mathbb{R}^{n_{x}}\right)$ a
state function and $u\in
L^{\infty}\left(\interval[openright]{0}{\infty},U\right)$ a control function,
consider the following continuous optimal control problem (OCP) on an infinite
time horizon:

Problem 1 (Continuous OCP).

Minimize

$J(x_{0},u):=\int_{0}^{\infty}{\exp(-\lambda t)C(x(t),u(t))dt}$ (1)

subject to the constraints

$\begin{array}[]{rcl}x^{\prime}(t)&=&F(x(t),u(t)),\\\ x(0)&=&x_{0},\\\
u(t)&\in&U.\end{array}$

Herein, $x_{0}$ is a given initial state, $F:\mathbb{R}^{n_{x}}\times
U\to\mathbb{R}^{n_{x}}$ is a continuous function describing the state
dynamics, and $C:\mathbb{R}^{n_{x}}\times U\to\mathbb{R}$ is a continuous
running cost function.

Let the value function $V$ of the optimal control problem (1) be defined as

$V(x_{0}):=\inf\\{J(x_{0},u)\mid u(t)\in U\ \text{for almost every}\
t\in\mathbb{R}_{\geq 0}\\}.$

It can be shown [12] that $V$ solves the following equation

$\lambda V(x)+\sup_{u\in U}\\{-\nabla V(x)^{T}F(x,u)-C(x,u)\\}=0,$

which is called the Hamilton-Jacobi-Bellman (HJB) equation [13].

For a given $x\in\mathbb{R}^{n_{x}}$, the feedback law can be then computed as

$u(x)=\underset{u\in U}{\arg\min}\\{\nabla V(x)^{T}F(x,u)+C(x,u)\\},$

whenever this minimum exists. Note that the gradient of $V$ is denoted by
$\nabla V(x)$ for $V$ supposed to be sufficiently smooth.

###  2.2 Discrete Time

Consider the following discrete optimal control problem (OCP) with a constant
discount rate $0<\beta\leq 1$:

Problem 2 (Discrete OCP).

Minimize

$j(x_{0},\\{u_{n}\\}_{n\in\mathbb{N}_{0}}):=\sum_{n=0}^{\infty}{\beta^{n}c(x_{n},u_{n})}$
(2)

subject to the constraints

$\begin{array}[]{ll}\begin{array}[]{rcl}x_{n+1}&=&f(x_{n},u_{n}),\\\
u_{n}&\in&U\end{array}&\text{for}\ n=0,1,\dots\end{array}$

Herein, $f:\mathbb{R}^{n_{x}}\times U\to\mathbb{R}^{n_{x}}$ is a continuous
function describing the discrete state dynamics, and
$c:\mathbb{R}^{n_{x}}\times U\to\mathbb{R}$ is a continuous running cost
function.

A detailed look at the objective function of problem (2) reveals that it is
constructed in such a way that the closer the summands to the initial state
$x_{0}$ are, the stronger their influence on $j$ is due to the weighting
parameters $\beta^{n},\ n\in\mathbb{N}_{0}$.

The value function $V_{h}$ of the discrete optimal control problem (2) is
defined as

$V_{h}(x_{0}):=\inf\\{j(x_{0},\\{u_{n}\\}_{n\in\mathbb{N}_{0}})\mid u_{n}\in
U\ \text{for}\ n\in\mathbb{N}_{0}\\}$ (3)

and represents the greatest lower bound of the total cost for all
trajectories, whose starting points are $x_{0}$ [14]. It can be shown [13]
that $V_{h}$ solves the following fixed-point equation for $V_{h}(x)$

$V_{h}(x)=\inf_{u\in U}\\{c(x,u)+\beta V_{h}(f(x,u))\\},$ (4)

which is a discrete version of the HJB equation. The derivation of the
equation exploits the definition of $V_{h}$ and involves, inter alia, the
insertion of equation (2) into equation (3). For reasons of space, it is
omitted here. The interested reader is instead referred to [14] for more
details.

Equation (4) can be solved using the following fixed-point iteration for an
initial guess $V_{h}^{0}$

$V_{h}^{m+1}(x)=\inf_{u\in U}\\{c(x,u)+\beta
V_{h}^{m}(f(x,u))\\}\quad\text{for}\ m\in\mathbb{N}_{0}$ (5)

until the sequence of functions $\\{V_{h}^{m}\\}_{m\in\mathbb{N}_{0}}$
converges [12].

For $\beta\in\interval[open]{0}{1}$, the operator

$T[V](x):=\inf_{u\in U}\\{c(x,u)+\beta V(f(x,u))\\}$

is a $\beta$-contraction. Thus, equation (5) possesses a unique fixed point
according to the Banach fixed-point theorem. However, it cannot be guaranteed
that the discrete HJB equation has a unique fixed point for $\beta=1$. An
elegant way of handling this issue is the so-called Kruz̆kov transform [15]
discussed later in Section 4.1.

######  Remark 1.

Note that the explicit dependency on time $t$ of a non-autonomous OCP can be
eliminated by introducing a new state variable

$\tau(t):=t.\\\ $

In the continuous case, the new state vector corresponds then to

$\tilde{x}(t)=\begin{bmatrix}x(t)\\\ \tau(t)\end{bmatrix},$

which solves the following initial value problem (IVP):

$\displaystyle\tilde{x}^{\prime}(t)$
$\displaystyle=\begin{bmatrix}F(\tau(t),x(t),u(t))\\\ 1\end{bmatrix},$
$\displaystyle\tilde{x}(0)$ $\displaystyle=\begin{bmatrix}x_{0}\\\
0\end{bmatrix}.$

Similarly, in the discrete case it holds

$\displaystyle\tilde{x}_{n+1}=\begin{bmatrix}x_{n+1}\\\
\tau_{n+1}\end{bmatrix}=\begin{bmatrix}f(x_{n},u_{n})\\\
\tau_{n}+h\end{bmatrix}\quad\text{for}\ n\in\mathbb{N}_{0},$

where $h>0$ is the step size of the time domain grid assumed to be
equidistant. Therefore, a value function $V_{h}(t,x(t))$ explicitly depending
on $t$ can be transformed into a function $V_{h}(\tilde{x}(t))$ depending
solely on the state variables.

###  2.3 Relation between Continuous and Discrete Optimal Control Problems

In order to show a relation between the continuous OCP and its discrete
counterpart, the continuous problem is discretized on the grid $t_{n}=nh,\
n\in\mathbb{N}_{0},$ with a step size $h>0,\ x_{n}:=x(t_{n})$ and
$u_{n}:=u(t_{n})$. This implies:

$J(x_{0},u)=h\sum_{n=0}^{\infty}{\exp(-\lambda
t_{n})C(x_{n},u_{n})}=\sum_{n=0}^{\infty}{{\exp(-\lambda
h)}^{n}hC(x_{n},u_{n})}.$

Hence, the relation holds with $\beta=\exp(-\lambda h)$ and $c(x,u)=hC(x,u)$.

To obtain a relation between the functions $f$ and $F$, the Euler method is
exploited for dealing with the ODE describing the state dynamics. This results
in

$F(x_{n},u_{n})=x^{\prime}(t_{n})\approx\frac{x_{n+1}-x_{n}}{h}=\frac{f(x_{n},u_{n})-x_{n}}{h}.$

Thus, it can be recognized that $f(x,u)=x+hF(x,u)$.

##  3 Problem Formulation

This section is devoted to the exact formulation of the problem. The overall
goal is to determine an optimal schedule for a small fleet of autonomous VTOLs
and to find the respective trajectories. Hereby, their dynamics and the
condition that while one VTOL executes its operation, all others have to idle
are taken into account. In our case, the mission for all VTOLs is to reach a
static target from some starting points while avoiding moving obstacles such
that the total mission time, i.e. the sum of the respective flight durations,
is minimal. Note that flight durations can be penalized.

This setting can be formulated as a mixed-integer bilevel optimization
problem, where the upper level problem deals with the scheduling of autonomous
VTOLs, while the lower level one delivers optimal trajectories. In what
follows, some notation used in this paper is introduced. Afterwards, the
underlying problem’s mathematical formulation is discussed in detail.

Throughout this work, the number of autonomous VTOLs is denoted by
$N\in\mathbb{N}$ with the corresponding index set
$\mathcal{N}:=\\{1,\dots,N\\}$. Moreover, $P\in\mathbb{N}_{0}$ and
$\mathcal{P}:=\\{1,\dots,P\\}$ indicate the number of (moving) obstacles,
whose dynamics is supposed to be known, and its index set, respectively. In
addition, the fixed starting positions of all aerial vehicles are assumed to
be given. The overall goal is to schedule VTOLs and to simultaneously compute
their optimal trajectories in some time-dependent domain
$\Omega(t)\subset\mathbb{R}^{n_{x}}$ with $t\in\mathbb{R}_{\geq 0}$ and
$n\in\\{2,3\\}$ such that all of them reach a predefined static target set
$\Omega_{T}\subset\Omega(t),\ t\in\mathbb{R}_{\geq 0},$ while avoiding
possible collisions with any of the obstacles. Note that the latter form the
time-dependent set $\Omega_{O}(t)$ with
$\Omega_{O}(t)=\underset{i\in\mathcal{P}}{\cup}\Omega_{O,i}(t)\subset\Omega(t),\
t\in\mathbb{R}_{\geq 0}$. The set, across which VTOLs are allowed to move at
any time point, is denoted by $\Omega_{S}(t)$. The whole domain $\Omega(t)$
can be thus represented as
$\Omega(t)=\Omega_{T}\cup\Omega_{O}(t)\cup\Omega_{S}(t),\ t\in\mathbb{R}_{\geq
0}$. An illustrative example of $\Omega(t_{\star})$ with two obstacles is
given in Figure 1 for a specific time point $t_{\star}\in\mathbb{R}_{\geq 0}$.

Figure 1: Example of $\Omega(t_{\star}),\ t_{\star}\in\mathbb{R}_{\geq 0}$. In
green: static target set $\Omega_{T}$, in orange: dynamic obstacles with
$\Omega_{O}(t_{\star})=\Omega_{O,1}(t_{\star})\cup\Omega_{O,2}(t_{\star})$, in
blue: time-dependent set $\Omega_{S}(t_{\star})$, across which VTOLs can move.

In addition, the following assumptions must apply within the framework of this
article. We briefly comment each of them.

######  Assumption 3.1.

The target region $\Omega_{T}$ can be reached by every VTOL in finite time.

This assumption is introduced in order to exclude situations where, for
instance, the target set is unreachable because of insurmountable obstacles.

######  Assumption 3.2.

The velocity of all VTOLs is assumed to be constant.

This assumption is legitimate, since a VTOL moves with a constant velocity
most of the flight. It allows to consider a simplified version of the
underlying OCPs, see Problem (7).

######  Assumption 3.3.

All obstacles stop their motion at some time point $t^{\star}$ bounded from
below by the time point, at which the last VTOL in the sequence reaches
$\Omega_{T}$.

This is a reasonable assumption, since the overall mission is said to be
fulfilled after the last aerial vehicle has reached the target set. It is
essential for our computational approach, see Section 4.1.

We are now at a position to formulate the whole problem with a bilevel
structure. As already mentioned before, the VTOLs in our consideration are
allowed to execute their respective operations only one by one. Hence, the
scheduling part of the problem can be formulated as the so-called job-shop
problem (JSP) [16].

For reasons of comprehensiveness, the bilevel problem with a coupling between
the two levels is presented as an upper and lower level problem and can be
formulated as

Problem 3 (Upper Level: Scheduling of VTOLs).

$\min_{t_{i},w_{ij}}\sum_{i=1}^{N}\left[t_{i}+\alpha D_{i}(t_{i})\right]$ (6a)
subject to $\displaystyle t_{i}+D_{i}(t_{i})-t_{j}\quad$
$\displaystyle\leq\quad$ $\displaystyle M(1-w_{ij}),$ (6b) $\displaystyle
t_{j}+D_{j}(t_{j})-t_{i}\quad$ $\displaystyle\leq\quad$ $\displaystyle
Mw_{ij},$ (6c) $\displaystyle t_{i}\quad$ $\displaystyle\in\quad$
$\displaystyle[t_{i,min},t_{i,max}],$ (6d) $\displaystyle w_{ij}\quad$
$\displaystyle\in\quad$ $\displaystyle\\{0,1\\},\ i,j\in\mathcal{N},i<j.$ (6e)

Herein, $t_{i}$ denotes the starting time of VTOL $i,\ i\in\mathcal{N},$
$\alpha$ is a flight duration penalization term, $M$ is some sufficiently
large constant, which guarantees that only one of constraints (6b), (6c) is
active at the same time. The binary variables $w_{i,j},\ i,j\in\mathcal{N},$
have the following meaning:

$w_{ij}=\begin{cases}1&\textnormal{if VTOL $i$ starts before VTOL $j$},\\\
0&\textnormal{otherwise}.\end{cases}$

Moreover, $D_{i}(t_{i}),\ i\in\mathcal{N},$ is a nonlinear function and
denotes the flight duration of the corresponding VTOL if it starts at time
$t_{i}$, which in turn is a continuous variable bounded by $t_{i,min}$ and
$t_{i,max}$ from below and above, respectively. Furthermore, the following
relation holds due to Assumption 3.2

$D_{i}(t_{i})=d(t_{i})/v_{i},\quad i\in\mathcal{N},$

where $v_{i}$ is the average velocity of VTOL $i$. For simplicity reasons, we
assume that $v_{i}=v,\ \forall i\in\mathcal{N}$. The numerator in the previous
equality is a parametric function, which denotes the length of the shortest
path to the target at time $t_{i},\ i\in\mathcal{N}$, and can be computed from
the following parametric optimal control problem:

Problem 4 (Lower Level: Shortest Path).

$d(t_{i}):=\min_{x_{i},u_{i}}\int_{t_{i}}^{\infty}{\chi_{\Omega_{S}(t)}(x_{i}(t))dt}$
(7a) subject to $\displaystyle x_{i}^{\prime}(t)$ $\displaystyle=\ $
$\displaystyle u_{i}(t)\quad$ $\displaystyle\textit{for a.e.}\
t\in\interval[openright]{t_{i}}{\infty},$ (7b) $\displaystyle x_{i}(t)$
$\displaystyle\in\ $ $\displaystyle\Omega_{S}(t)\cup\Omega_{T}\quad$
$\displaystyle\textit{for every}\ t\in\interval[openright]{t_{i}}{\infty},$
(7c) $\displaystyle u_{i}(t)$ $\displaystyle\in\ $ $\displaystyle
U_{i}(x_{i})=\begin{cases}\\{u\in\mathbb{R}^{n_{u}}\mid||u||=1\\},\
x_{i}\in\Omega_{S}(t)\\\ \\{0\in\mathbb{R}^{n_{u}}\\},\
x_{i}\in\Omega_{T}\end{cases}\quad$ $\displaystyle\textit{for a.e.}\
t\in\interval[openright]{t_{i}}{\infty},$ (7d) $\displaystyle x_{i}(t_{i})$
$\displaystyle=\ $ $\displaystyle\overline{x}_{i},\quad$ (7e)

where $\overline{x}_{i}$ is the initial position of VTOL $i$, and $\chi_{A}$
is the indicator function defined as

$\chi_{A}(x)=\begin{cases}1&\text{if}\ x\in A,\\\
0&\text{otherwise}.\end{cases}$

Note that in case that $n_{u}=2$, the control sets $U_{i}(x_{i})$ for
$x_{i}\in\Omega_{S}(t),\ i\in\mathcal{N},$ can be represented in the following
way using polar coordinates:

$U_{i}(x_{i})=\left\\{\left(\ \begin{matrix}u_{1}\\\ u_{2}\end{matrix}\
\right)\ \Big{|}\ u_{1}=\cos{\alpha_{i}},\ u_{2}=\sin{\alpha_{i}},\
\alpha_{i}\in\interval[openright]{0}{2\pi}\right\\}.$

Similarly, the spherical coordinate system can be exploited for describing
$U_{i}(x_{i})$ for $x_{i}\in\Omega_{S}(t),\ i\in\mathcal{N},$ if $n_{u}=3$:

$U_{i}(x_{i})=\left\\{\left(\ \begin{matrix}u_{1}\\\ u_{2}\\\
u_{3}\end{matrix}\ \right)\ \bigg{|}\ u_{1}=\cos{\alpha_{i}}\sin{\beta_{i}},\
u_{2}=\sin{\alpha_{i}}\sin{\beta_{i}},\
u_{3}=\cos{\beta_{i}}\,\alpha_{i}\in\interval[openright]{0}{2\pi},\
\beta_{i}\in\interval{0}{\pi}\right\\}.$

By setting

$\displaystyle F(x,u)$ $\displaystyle=u,$ $\displaystyle C(x,u)$
$\displaystyle=\chi_{\Omega_{S}(t)}(x),$ $\displaystyle\lambda$
$\displaystyle=0,$

it becomes clear that this problem fits into the framework of OCP classes, see
(2.1). Hence, it can be solved with the value function approach. This is very
beneficial in our case, since only one value function has to be computed for
all problems due to their structure. In this case, it has to be computed for
$t\in\mathbb{R}_{t\geq t_{min}},\ t_{min}=\min_{i\in\mathcal{N}}t_{i,min}$.
Without loss of generality, we will hereinafter always assume $t_{min}=0$.

Obviously, problem (7) is time-dependent because the underlying domain depends
on time. In order to fit this problem into the framework of OCP classes, steps
from Remark 1 have to be exploited. However, they are omitted here due to
readability reasons and appear only in the discrete case, see Section 4.1. It
can be shown that solving problem (7) is equivalent to solving the Eikonal
equation given by:

$\begin{split}||\nabla V_{h}(x_{i})||&=1\quad\text{for}\
x_{i}\in\Omega_{S}(t),\ t\in\mathbb{R}_{t\geq 0},\\\
V_{h}(x_{i})&=0\quad\text{for}\ x\in\partial\Omega_{T},t\in\mathbb{R}_{t\geq
0}.\end{split}$

Therefore, the solution of this problem yields the shortest paths from
$\overline{x}_{i},\ i\in\mathcal{N},$ to the target set $\Omega_{T}$.

##  4 Solution Procedure

This section highlights and explains main steps required to solve the bilevel
mixed-integer optimization problem (6) - (7). The solution procedure can be
divided into three main steps indicated below:

  1. 1.

Solve the lower level problems corresponding to the parametric optimal control
problems (7). Its output, as discussed above, is a single value function
$V_{h}(t,x)$ depending on both the starting time $t$ and the starting position
$x$.

  2. 2.

Insert all given VTOLs’ starting positions $\overline{x}_{i},\
i\in\mathcal{N},$ into the value function $V_{h}$ obtained in step (1) and set
$d_{i}(t_{i}):=V_{h}(t_{i},\overline{x}_{i})$ to obtain time-dependent
shortest-path functions. The bilevel problem (6) - (7) is hence recast into a
single-level MINLP. It is then piecewise linearized and solved by a robust and
efficient MILP-solver. A solution to that problem is an optimal sequence of
VTOLs provided by the optimal starting times $t^{\ast}_{i},\ i\in\mathcal{N}$.

  3. 3.

Compute the optimal trajectory of every VTOL via the forward recursion.

###  4.1 First Step: Solving the Optimal Control Problem

In the first step of the solution procedure, the parametric optimal control
problems (7) need to be solved. Before we go deeper into detail, let us first
explain the idea behind the formulation of (7). Starting at a point
$\overline{x}_{i}$, a circle with a radius of one and a midpoint at
$\overline{x}_{i}$ is considered. The aim is to find that point lying on the
boundary of this circle, which has the minimum distance to the target set,
avoiding a collision with dynamic obstacles. The distance between
$\overline{x}_{i}$ and the newly found point is thus equal to one. In the next
step, a new circle with its midpoint lying in this new point is considered.
The whole procedure is repeated until the last circle crosses the target area.
In the end, the length of the shortest path from point $\overline{x}_{i}$ to
the target set $\Omega_{T}$ can be computed as a sum of ones. Figure 2
visualizes the first steps during the computation procedure according to this
idea.

Figure 2: Approach for the shortest path computation. In green: static target
set $\Omega_{T}$, in orange: a dynamic obstacle $\Omega_{O}(t)$, in purple: a
moving VTOL

To obtain a solution of problem (7), the value function approach discussed in
Section 2.2 is exploited. Although it is known to suffer from the "curse of
dimensions", in the case of the shortest-path problem the dimensionality is
characterized only by the time and space variables. Thus, the problem
possesses minimum possible dimensionality. Moreover, since the value function
is only required to obtain the optimal starting times, simple VTOL models can
be used to determine them. Afterwards, more complex models can be exploited in
order to compute more precise trajectories. Note that the corresponding OCPs
will be no longer parametric in this case and can be therefore solved with
standard techniques.

To determine the setting of the discrete problem, we use the relation between
it and its continuous counterpart explained in Section 2.3 and the argument
from Remark 1. As a result, we obtain:

$\begin{array}[]{rcl}f(x,u)&=&x+hu,\\\ c(x,u)&=&\begin{cases}h&\text{if}\
x\in\Omega_{S}(t_{k}),\\\ 0&\text{otherwise},\end{cases}\\\ \beta&=&1,\\\
t_{k+1}&=&t_{k}+h,\quad k\in\mathbb{N}_{0}.\end{array}$

For the value function, it holds

$V_{h}(t_{k},x)=\begin{cases}\infty&\text{for}\
x\in\left(\mathbb{R}^{n_{x}}\setminus\Omega\right)\cup\Omega_{O}(t_{k}),\\\
0&\text{for}\ x\in\Omega_{T}\end{cases}\quad\text{for}\ k\in\mathbb{N}_{0}.$
(8)

The value function can be computed iteratively using formula (5) indicated in
Section 2. However, this might cause several difficulties. From a theoretical
point of view, it is not clear a priori whether equation (4) possesses a
unique solution for $\beta=1$. From a numerical point of view, the fact that
$V_{h}$ can attain infinite values poses another problem, since they might be
propagated [14]. An effective remedy is provided by the so-called Kruz̆kov
transform [15, 14] defined as:

$v_{h}(t_{k},x):=\exp(-V_{h}(t_{k},x)),\quad k\in\mathbb{N}_{0}.$ (9)

Applying (9) to (4) and taking into account (8) yields

$v_{h}(t_{k},x)=\begin{cases}0&\text{for}\
x\in\left(\mathbb{R}^{n_{x}}\setminus\Omega\right)\cup\Omega_{O}(t_{k}),\\\
\underset{u\in U}{\sup}\\{\exp(-c(x,u))v_{h}(t_{k+1},f(x,u))\\}&\text{for}\
x\in\Omega_{S}(t_{k}),\\\ 1&\text{for}\ x\in\Omega_{T}\end{cases}$

for $k\in\mathbb{N}_{0}$.

Similarly to (5), the transformed value function can be computed iteratively

$v_{h}^{m+1}(t_{k},x):=\sup_{u\in
U}\\{\exp(-c(x,u))v_{h}^{m}(t_{k+1},f(x,u))\\}\quad\text{for}\
k,m\in\mathbb{N}_{0}$

with the following initial guess

$v_{h}^{0}(t_{k},x)=\begin{cases}1&\text{if}\ x\in\Omega_{T},\\\
0&\text{otherwise}\end{cases}\quad\text{for}\ k\in\mathbb{N}_{0}.$

Provided that $c(x,u)$ is positive, the sequence of functions
$\\{v_{h}^{m}\\}_{m\in\mathbb{N}_{0}}$ will converge to some function
$\hat{v}_{h}$. The original value function can be then obtained using the
following backtransformation:

$V_{h}(t_{k},x)=-\ln\big{(}\hat{v}_{h}(t_{k},x)\big{)},\quad
k\in\mathbb{N}_{0}.$

One of problems, which arises during the numerical computation of the value
function $V_{h}$, relates to the question how the infinite dimensional
problems can be handled in practice. To overcome this issue, the following two
cases are considered:

  1. 1.

All obstacles are supposed to move periodically with some global period $T$.
In this case, the value function can be computed in the following way:

$\displaystyle v_{h}^{m+1}(t_{k},x)$ $\displaystyle=\sup_{u\in
U}\\{\exp(-c(x,u))v_{h}^{m}(t_{k+1},f(x,u))\\}$ $\displaystyle\text{for}\
k=0,\dots,K-1,$ $\displaystyle v_{h}^{m+1}(t_{K},x)$ $\displaystyle=\sup_{u\in
U}\\{\exp(-c(x,u))v_{h}^{m}(t_{0},f(x,u))\\},$

where $Kh=T$.

  2. 2.

Due to Assumption 3.3, the underlying domain can be viewed as static after
some time point $t^{\star}$, hence

$\displaystyle v_{h}^{m+1}(t_{k},x)$ $\displaystyle=\sup_{u\in
U}\\{\exp(-c(x,u))v_{h}^{m}(t_{k+1},f(x,u))\\}$ $\displaystyle\text{for}\
k=0,\dots,K-1,$ $\displaystyle v_{h}^{m+1}(t_{K},x)$ $\displaystyle=\sup_{u\in
U}\\{\exp(-c(x,u))v_{h}^{m}(t_{K},f(x,u))\\},$

where $Kh>t^{\star}$.

In numerical computations, the space domain is discretized with some
equidistant gridpoints. To obtain values of the value function in non-grid
points, piecewise linearization can be exploited.

###  4.2 Second Step: Solving the Scheduling Problem

Using the value function obtained in the previous step, the duration functions
of each VTOL can be now defined as

$d_{i}(t_{i}):=V_{h}(t_{i},\overline{x}_{i}),\quad i\in\mathcal{N}.$

Recall that $\overline{x}_{i}$ is the initial starting position of VTOL $i,\
i\in\mathcal{N}$. Hence, one ends up with MINLP (6).
In general, MINLPs are very hard to solve, especially for very nonlinear
problems. To overcome this issue, this work proposes an approximation
technique with piecewise linear functions. It eventually yields an MILP and
allows to track the linearzation error.

To this end, define a closed interval in $\mathbb{R}$ by

$I_{a}^{b}:=\interval{a}{b}.$

For time intervals $I_{t_{i,min}}^{t_{i,max}},\ i\in\mathcal{N},$ consider
their discretizations

$\mathcal{G}_{i}:=\\{t_{i,1},t_{i,2},\dots,t_{i,K_{i}}\\},$

where $t_{i,1}=t_{i,min},\ t_{i,K_{i}}=t_{i,max}$, and $K_{i}$ denotes the
number of discretization points in the grid $\mathcal{G}_{i},\
i\in\mathcal{N}$.

Next, consider a nonlinear function $g:\interval{t_{a}}{t_{b}}\to\mathbb{R}$
and its linear approximation $\psi(t)$ satisfying the conditions
$g(t_{a})=\psi(t_{a})$ and $g(t_{b})=\psi(t_{b})$. The maximum underestimator
error is then defined as

$e_{u}(g,I_{t_{a}}^{t_{b}}):=\max_{t\in I_{t_{a}}^{t_{b}}}g(t)-\psi(t),$

whereas the maximum overestimator error is indicated by

$e_{o}(g,I_{t_{a}}^{t_{b}}):=\max_{t\in I_{t_{a}}^{t_{b}}}\psi(t)-g(t).$

For more details, the interested reader is referred to [8] and [17].

Figure 3 shows a piecewise linearization of a nonlinear function $g(t_{i})$,
whereby the so-called envelopes around the piecewise linear function are
highlighted in orange.

Figure 3: Piecewise linearization with envelope functions

In order to obtain a fully linear problem, which can be directly handled by a
linear solver, two sets of auxiliary variables are additionally introduced.
Namely, binary variables $z_{i,k}$ and real variables $\lambda_{i,k}$ for
$i\in\\{1,\dots,N\\}$ with $k\in\\{1,\dots,K_{i}-1\\}$ and
$k\in\\{1,\dots,K_{i}\\}$, respectively. The integer variables $z_{i,k}$ are
called pointers and indicate that subinterval, within which the optimal value
of $t_{i}$ lies. As to the real variables $\lambda_{i,k}$, they represent the
coefficients in the optimal subinterval required to build a convex combination
there. Obviously, $t_{i}$ attains its optimal value on exactly one
subinterval, wherefore there is only one non-zero pointer.

These arguments lead to the following constraints:

$\displaystyle\begin{split}&\sum_{k=1}^{K_{i}-1}z_{i,k}=1,\\\
&z_{i,k}\in\left\\{0,1\right\\},\quad k\in\\{1,\dots,K_{i}-1\\}\end{split}$
(10)

and

$\displaystyle\begin{split}\sum_{k=1}^{K_{i}}\lambda_{i,k}=1.\end{split}$ (11)

In addition, the following inequalities have to be incorporated into the model
to guarantee that two consecutive lambdas are unequal to zero:

$\displaystyle\begin{split}\lambda_{i,1}&\leq z_{i,1},\\\ \lambda_{i,k}&\leq
z_{i,k-1}+z_{i,k},\quad k\in\\{2,\dots,K_{i}-1\\},\\\ \lambda_{i,K_{i}}&\leq
z_{i,K_{i}-1}.\\\ \end{split}$ (12)

Thus, the following additional constraints result from the preceding
considerations regarding the piecewise linearization technique for
$i\in\\{1,\dots,N\\}$:

$\displaystyle t_{i}$ $\displaystyle=\sum_{k=1}^{K_{i}}\ \lambda_{i,k}\
t_{i,k},$ (13a) $\displaystyle D_{i,pw}$ $\displaystyle:=\sum_{k=1}^{K_{i}}\
\lambda_{i,k}\ D_{i}(t_{i,k})+e_{i},$ (13b) $\displaystyle e_{i}$
$\displaystyle\leq\sum_{k=1}^{K_{i}-1}e_{u}\left(D_{i},I_{t_{i,k}}^{t_{i,k+1}}\right)z_{i,k},$
(13c) $\displaystyle e_{i}$
$\displaystyle\geq-\sum_{k=1}^{K_{i}-1}e_{o}\left(D_{i},I_{t_{i,k}}^{t_{i,k+1}}\right)z_{i,k}.$
(13d)

Herein, (13a) represents the starting time variable $t_{i}$ as a linear
combination of discretization points, (13b) corresponds to the linearizations
of the nonlinearities $D_{i}(\cdot)$ with the corresponding envelope terms.
Those are bounded from above and below by the maximum under- and
overestimators (13c)-(13d), respectively.

With this in mind, the mixed-integer linear relaxation of problem (6) can be
finally formulated as

Problem 5 (MILP).

$\begin{array}[]{l}\displaystyle\min_{t_{i},w_{ij},e_{i}}\
\sum_{i=1}^{N}\left[t_{i}+\alpha D_{i,pw}\right],\\\
\textnormal{s.t.}\quad\eqref{constr:sch1}-\eqref{constr:wij},\eqref{eq:PW_contraints_z}-\eqref{eq:PW_all}.\end{array}$
($\Pi$)

Herein, all continuous variables in constraints (6b) - (6d) are replaced by
their discrete counterparts.

The formulation of the adaptive numerical algorithm used for the refinement of
MILP relaxations is inspired by similar algorithms from [8] and [18]. The idea
is outlined in Algorithm 1 and can be briefly described in the following way.
First, some coarse discretizations $\mathcal{G}_{i}$ of the corresponding time
intervals $I_{t_{i,min}}^{t_{i,max}},\ i\in\mathcal{N},$ are chosen. The
algorithm then detects the subintervals where a predefined error bound is
violated with the aim of refining them. In the next iteration step, the new
refined time discretizations are examined. The whole procedure is repeated
until the error bound is satisfied for all $i\in\mathcal{N}$ [8]. As to the
termination condition, it was proven in [18] that the number of iteration
steps is guaranteed to be finite provided that problem (6) is feasible, and
all nonlinear functions are continuous.

Algorithm 1 Global optimization through adaptive refinement of MILP
relaxations

1:An MINLP (6) , initial tolerance $\epsilon_{0}$, the maximal linearization
error $\epsilon$ and the number of new discretization points $L$.

2:If (6) is feasible, the algorithm returns an optimal solution
$t^{\ast}_{i},\ i\in\mathcal{N},$ of ($\Pi$), which is an MILP relaxation of
(6), such that $\lvert D_{i}(t^{\ast}_{i})-D_{i,pw}\rvert\leq\varepsilon$,
$i\in\mathcal{N},$ and the corresponding objective value is minimal for any
admissible point of (6).

3:Choose initial discretizations $\mathcal{G}^{0}_{i},\ i\in\mathcal{N},$ such
that the corresponding nonlinearities do not exceed the initial tolerance
value.

4:Set $r\leftarrow 0$

5:repeat

6: Construct an MILP relaxation $\Pi^{r}$ of (6) from $\mathcal{G}^{r}_{i},\
i\in\mathcal{N}$.

7: Solve $\Pi^{r}$.

8: if $\Pi^{r}$ is feasible then

9: Set $t_{i}^{r},\ i\in\mathcal{N}\leftarrow$ optimal solution of $\Pi^{r}$

10: Set $D^{r}_{i,pw},\ i\in\mathcal{N}\leftarrow$ linear approximation of
nonlinear function $D_{i}(t_{i}^{r})$
---

11: else

12: return infeasible

13: end if

14: Set stop $\leftarrow$ true

15: for all $i\in\mathcal{N}$ do

16: Set $err^{r}\leftarrow\lvert D_{i}(t^{r}_{i})-D^{r}_{i,pw}\rvert$

17: if $err^{r}>\epsilon$ then

18: Set $n_{0}\leftarrow$ the index of the corresponding nonzero pointer
$z_{i,n_{0}}^{r}$ from the solution of $\Pi^{r}$
---

19: Set $h^{r}\leftarrow\frac{t^{r}_{i,n_{0}+1}-t^{r}_{i,n_{0}}}{L+1}$
---

20: Set
$\mathcal{G}^{r}_{new}\leftarrow\\{t^{r}_{i,n_{0}}+h^{r},t^{r}_{i,n_{0}}+2h^{r},\dots,t^{r}_{i,n_{0}}+Lh^{r}\\}$
---

21: Set
$\mathcal{G}_{i}^{r+1}\leftarrow\mathcal{G}_{i}^{r}\cup\mathcal{G}^{r}_{new}$
---

22: Set stop $\leftarrow$ false

23: else

24: Set $\mathcal{G}_{i}^{r+1}\leftarrow\mathcal{G}_{i}^{r}$

25: end if

26: end for

27: Set $r\leftarrow r+1$

28:until stop

29:return $t^{r-1}_{i},\ i\in\mathcal{N}$

###  4.3 Third Step: Computation of Optimal Trajectories

In the last step, the optimal trajectories are computed, see [13]. First, the
optimal starting times $t^{\ast}_{i},\ i\in\mathcal{N},$ are rounded up to the
next node in the time grid. This step is not necessary, but is performed here
due to safety reasons in order to avoid motions’ overlaps. Without performing
this step, interpolation of the value function in time is required in
addition. We denote those new starting times by $t_{k_{i}},\ i\in\mathcal{N}$.
For each $\ i\in\mathcal{N}$, optimal trajectories can be computed as follows:

$\displaystyle u^{\ast}_{i}(t_{k})$ $\displaystyle=\underset{u\in
U}{\arg\min}\\{h+V_{h}(t_{k+1},f(x^{\ast}_{i}(t_{k}),u))\\},$ $\displaystyle
x^{\ast}_{i}(t_{k+1})$
$\displaystyle=f(x^{\ast}_{i}(t_{k}),u^{\ast}_{i}(t_{k})),\quad k=k_{i},\
k_{i}+1,\ k_{i}+2,...\ .$

If $V_{h}(t_{k+1},f(x^{\ast}_{i}(t_{k}),u))=0$ for some $k$,
$x^{\ast}_{i}(t_{k+1})$ is the last point in the trajectory meaning that the
algorithm terminates.

Note that error estimates are possible, and that the error in $V_{h}$ can be
controlled.

##  5 Numerical Results

All numerical experiments were conducted on an Intel(R) Core i7 processor with
a 2.80GHz-CPU and an 8GB-RAM using the robust and efficient GUROBI MILP-solver
[19] based on the Branch-and-Bound algorithm. For reasons of space, the idea
behind this algorithm is not explained here. The interested reader is instead
referred to [20].

Our approach was extensively tested for up to $N=8$ VTOLs to show its
feasibility. A periodic problem was considered. The time-space domain $\Omega$
with $n=2$ and the control set $U$ were replaced by their finite
discretizations. The value function $V_{h}$ was then obtained by its pointwise
evaluation at the corresponding discretization points. All results presented
below relate to the problem setting described as follows.

In the state space
$\Omega=\interval{-200}{200}\times\interval{-200}{200}\subset\mathbb{R}^{2}$
with $\Omega_{T}$ lying in its center, we consider $N=8$ VTOLs and $P=4$
dynamic obstacles. VTOLs in our consideration have the following initial
positions:

$\displaystyle\overline{x}_{1}$
$\displaystyle=\overline{x}_{5}=\\{-192,-192\\}$\mathrm{m}$,$
$\displaystyle\overline{x}_{2}$
$\displaystyle=\overline{x}_{6}=\\{192,-192\\}$\mathrm{m}$,$
$\displaystyle\overline{x}_{3}$
$\displaystyle=\overline{x}_{7}=\\{192,192\\}$\mathrm{m}$,$
$\displaystyle\overline{x}_{4}$
$\displaystyle=\overline{x}_{8}=\\{-192,192\\}$\mathrm{m}$.$

They move with a constant velocity $v=10$\mathrm{m}$/$\mathrm{s}$$. The
circular obstacles, with $R_{O}=64$\mathrm{m}$$ radius each, move periodically
around the target region $\Omega_{T}$ along a circular trajectory with radius
$R_{C}=100$\mathrm{m}$$. The initial positions of their centers are given by
$\overline{x}_{1}=[R_{C},0]^{T},\overline{x}_{2}=[0,R_{C}]^{T},\overline{x}_{3}=[-R_{C},0]^{T}$
and $\overline{x}_{4}=[0,-R_{C}]^{T}$, respectively. Having a period of
$80$\mathrm{s}$$ (full rotation around the target), their velocity, which is
assumed to be constant for convenience, is approximately equal to
$7.85$\mathrm{m}$/$\mathrm{s}$$.

Note that the value function $V_{h}(t,x)$ evaluated at some point
$\overline{x}$ corresponds to a time-dependent function, which indicates the
length of the shortest path from that point to the target region $\Omega_{T}$.
Figure 4 illustrates this function for a VTOL starting at
$\overline{x}=[192,192]^{T}$\mathrm{m}$$. Since all four obstacles move
periodically with a constant and identical velocity, it suffices to analyze
the plot during the first 20$\mathrm{s}$. The shortest path is observed when
the VTOL starts moving after $2$\mathrm{s}$$, while the longest path has a
length of nearly 315$\mathrm{m}$ and refers to the case where VTOL commences
its motion at time $t=11$\mathrm{s}$$. This observation can be explained by
the fact that the VTOL does not collide with any of the obstacles if it starts
at $t=2$\mathrm{s}$$ and moves straight towards the target. On the contrary,
if the VTOL starts at $t=11$\mathrm{s}$$, it has to avoid a collision with an
obstacle and therefore flies around it. For that reason, the path is in the
second case much longer than in the first one. After 20$\mathrm{s}$, the next
obstacle takes place of the previous one, and the same arguments apply.

Figure 4: Shortest-path function of a VTOL starting at point
$\overline{x}=[192,192]^{T}$\mathrm{m}$$

Next, we present a color map shown in Figure 5, which visualizes the value
function evaluated at the initial time point $t_{0}=0$. To be more specific,
for any point
$x\in\Omega(t_{0})=\interval{-200}{200}\times\interval{-200}{200}$ this plot
indicates the length of the shortest path (taking into account obstacle
avoidance) from a VTOL located at this point to the target region $\Omega_{T}$
with radius $R_{T}=10$\mathrm{m}$$. Obviously, if the initial position of a
VTOL lies somewhere within any of the obstacles, the value function $V_{h}$ is
equal to an infinite value according to our convention. Therefore, the four
circles in Figure 5 are colored yellow. Similarly, if the VTOL starts at a
point lying within the target set, its distance to this set is zero. This
explains why $\Omega_{T}$ illustrated as a circle in the middle of the color
map has the darkest color. The last interesting observation arises from the
semicircular regions depicted in Figure 5, which take such a shape because of
the obstacles’ round shape. If the initial position of a VTOL lies somewhere
within such blue region, the path from the VTOL to $\Omega_{T}$ would be
rather short due to the possibility of a collision-free flight along a nearly
straight trajectory.

Figure 5: Value function $V_{h}(t,x)$ evaluated at time $t_{0}=0$

In order to determine the average computational time required by the GUROBI
MILP solver to solve problem ($\Pi$), 10 experiments were run with the same
parameter setting. For $\alpha=1$, the solver needed on average
19.34$\mathrm{s}$ to find a numerical solution of the final MILP $\Pi$. In
case that $\alpha=10$, the computational time reduced drastically to
3.65$\mathrm{s}$. The achieved results yield that the solver runtime decreases
with a growing $\alpha$. This relation is caused by the fact that larger
values of the penalization term yield smaller search trees in the Branch-and-
Bound algorithm and hence shorter computational time.

Optimal schedules obtained by the solver for $\alpha=1$ and $\alpha=10$ are
shown in Tables 1 and 2, respectively. A closer look at the last column of
Tables 1 and 2 reveals that the larger the $\alpha$ values are, the smaller
the flight duration of a single VTOL becomes. This observation is totally
plausible, since the VTOLs are penalized more in the second case and therefore
try to complete their motions as quickly as possible. While the duration of a
single mission decreases with a growing $\alpha$, the total task duration
shows an opposite tendency. In our problem setting, the last VTOL reached the
target set after $225.57$\mathrm{s}$$ for $\alpha=1$, whereas the overall
mission was completed only after $247.885$\mathrm{s}$$ for $\alpha=10$. This
phenomenon can be explained by the fact that each VTOL generally tends to wait
longer when $\alpha$ increases in order to reduce its own flight duration. As
a result, the overall mission time increases.

Starting Number | Original Number | Starting Time [$\mathrm{s}$] | End Time [$\mathrm{s}$] | Duration [$\mathrm{s}$]
---|---|---|---|---
1 | 4 | 0 | 27.0094 | 27.0094
2 | 3 | 27.0094 | 55.1403 | 28.131
3 | 1 | 57.1793 | 84.9432 | 27.7639
4 | 2 | 84.9432 | 111.925 | 26.9819
5 | 5 | 111.925 | 142.105 | 30.1801
6 | 8 | 142.105 | 168.331 | 26.2255
7 | 6 | 168.331 | 197.768 | 29.4372
8 | 7 | 197.768 | 225.57 | 27.8021
Table 1: Optimal Schedule for $N=8$ VTOLs with $\alpha=1$ Starting Number | Original Number | Starting Time [$\mathrm{s}$] | End Time [$\mathrm{s}$] | Duration [$\mathrm{s}$]
---|---|---|---|---
1 | 1 | 2.63158 | 29.2895 | 26.658
2 | 3 | 39.4737 | 66.4896 | 27.0159
3 | 4 | 66.4896 | 94.339 | 27.8494
4 | 2 | 98.9736 | 126.316 | 27.3422
5 | 8 | 126.316 | 153.973 | 27.6573
6 | 6 | 158.427 | 185.965 | 27.5377
7 | 5 | 185.965 | 213.464 | 27.4995
8 | 7 | 221.053 | 247.855 | 26.8026
Table 2: Optimal Schedule for $N=8$ VTOLs with $\alpha=10$

Finally, a 2D visualization was created in order to illustrate the results
from above for two different values of the penalization term, namely for
$\alpha=1$ and $\alpha=10$. The respective videos can be found under:

https://youtu.be/QgMEM1QQP9M
https://youtu.be/BQ39fRMV-ew

##  6 Conclusion

The underlying work presented an efficient numerical framework for scheduling
of multiple autonomous VTOLs and their simultaneous optimal trajectory
planning. Formulated as a bilevel problem, the model is recast into a single-
level one by means of a value function approach. This technique involves such
advanced concepts as the Hamilton-Jacobi-Bellman equation originating from
dynamic programming theory and the Kruz̆kov transformation. The resulting
MINLP was piecewise linearized with the help of an elegant technique allowing
to track the linearization error. The solution of the final MILP was then
determined numerically by the GUROBI solver based on the Branch-and-Bound
algorithm. The numerical experiments yield the feasibility of the framework
presented in this article.

Possible model extensions encompass:

  * •

Extension of the optimal control problem at the lower level to cope with
multi-phase systems.

  * •

Consideration of alternative scheduling mixed-integer problems at the upper
level as, for example, vehicle routing.

  * •

Integration of other objectives at the upper level with the aim of maximizing
the mission success, minimizing the fuel consumption etc.

  * •

Integration of coupling effects by introducing common constraints and/or
common objectives.

  * •

Consideration of multi-stage problems with early starts.

## Funding

This material is based upon work supported, inter alia, by
dtec.bw–Digitalization and Technology Research Center of the Bundeswehr
[MissionLab] and by the Air Force Office of Scientific Research, Air Force
Materiel Command, USAF, under Award No. FA8655-20-1-7026. Any opinions,
findings, and conclusions or recommendations expressed in this publication are
those of the author(s) and do not necessarily reflect the views of the
funders.

## References

  * [1] EASA. Vertical take-off & landing [https://www.easa.europa.eu/en/domains/rotorcraft-vtol/VTOL]; n.d. Accessed: 2022-09-02.
  * [2] Kotoky A. Singapore’s tiny second airport eyes future as a flying taxi hub [https://www.bloomberg.com/news/articles/2022-07-04/singapore-s-tiny-second-airport-eyes-future-as-a-flying-taxi-hub]; 2022\. Accessed: 2022-09-08.
  * [3] Parzani C, Puechmorel S. On a hamilton-jacobi-bellman approach for coordinated optimal aircraft trajectories planning. Optimal Control Applications and Methods. 2017 12;39(2):933–948.
  * [4] Coupechoux M, Darbon J, Kélif JM, et al. Optimal trajectories of a uav base station using hamilton-jacobi equations ; 2021.
  * [5] Rogovs S, Nikitina V, Gerdts M. A novel mixed-integer programming approach for the aircraft landing problem. Frontiers in Future Transportation. 2022 09;3.
  * [6] Chandran BG, Balakrishnan H. A dynamic programming algorithm for robust runway scheduling. 2007 American Control Conference. 2007;:1161–1166.
  * [7] Mahdavi A, Carvalho M. Optimal trajectory and schedule planning for autonomous guided vehicles in flexible manufacturing system. In: 2018 Second IEEE International Conference on Robotic Computing (IRC); 01. IEEE; 2018.
  * [8] Gerdts M, Rogovs S, Valenti G. Solving complex intersection management problems using bi-level MINLPs and piecewise linearization techniques. Springer International Publishing; 2022. p. 255–273.
  * [9] Palagachev KD, Gerdts M. Numerical approaches towards bilevel optimal control problems with scheduling tasks. In: Math for the digital factory. Springer International Publishing; 2017. p. 205–228.
  * [10] Zanlongo SA, Abodo F, Long P, et al. Multi-robot scheduling and path-planning for non-overlapping operator attention. In: 2018 Second IEEE International Conference on Robotic Computing (IRC); 01. IEEE; 2018.
  * [11] Samá M, D'Ariano A, Pacciarelli D, et al. Optimal aircraft scheduling and flight trajectory in terminal control areas. In: 2017 5th IEEE International Conference on Models and Technologies for Intelligent Transportation Systems (MT-ITS); 06. IEEE; 2017.
  * [12] Bardi M, Capuzzo-Dolcetta I. Optimal control and viscosity solutions of hamilton-jacobi-bellman equations. Birkhäuser Boston; 1997.
  * [13] Bellman R. Dynamic programming. 1st ed. Princeton, NJ, USA: Princeton University Press; 1957.
  * [14] Schreiber AJE. Dynamic programming with radial basis functions and shepard’s method [dissertation]. Technische Universität München, Zentrum Mathematik, Wissenschaftliches Rechnen; 2016. PhD thesis.
  * [15] Kružkov SN. Generalized solutions of the hamilton-jacobi equations of eikonal type. i. formulation of the problems; existence, uniqueness and stability theorems; some properties of the solutions. Mathematics of the USSR-Sbornik. 1975 04;27(3):406.
  * [16] Błażewicz J, Domschke W, Pesch E. The job shop scheduling problem: Conventional and new solution techniques. European Journal of Operational Research. 1996;93(1):1–33.
  * [17] Burlacu R. Adaptive mixed-integer refinements for solving nonlinear problems with discrete decisions [dissertation]. Friedrich-Alexander-Universität Erlangen-Nürnberg; 2019. PhD thesis.
  * [18] Burlacu R, Geißler B, Schewe L. Solving mixed-integer nonlinear programmes using adaptively refined mixed-integer linear programmes. Optimization Methods and Software. 2019 01;35(1):37–64. Available from: https://doi.org/10.1080/10556788.2018.1556661.
  * [19] Gurobi Optimization LLC. Gurobi optimizer reference manual [https://www.gurobi.com/wp-content/plugins/hd_documentations/documentation/9.1/refman.pdf]; 2020\. Accessed: 2022-02-10.
  * [20] Winston W, Goldberg J. Operations research: Applications and algorithms. Belmont, CA: Thomson Brooks/Cole; 2004.