Karlsruher Schriften zur Anthropomatik

Band 58

Band 58

J. Pfrommer

Distr. Planning for Self-Organizing Production Systems

Julius Pfrommer

**Distributed Planning for Self-Organizing Production Systems**

Julius Pfrommer

## **Distributed Planning for Self-Organizing Production Systems**

Karlsruher Schriften zur Anthropomatik Band 58 Herausgeber: Prof. Dr.-Ing. habil. Jürgen Beyerer

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buchs.

# **Distributed Planning for Self-Organizing Production Systems**

by Julius Pfrommer

Karlsruher Schriften zur Anthropomatik

Herausgeber: Prof. Dr.-Ing. habil. Jürgen Beyerer

Eine Übersicht aller bisher in dieser Schriftenreihe erschienenen Bände finden Sie am Ende des Buchs.

Band 58

Karlsruher Institut für Technologie Institut für Anthropomatik und Robotik

Distributed Planning for Self-Organizing Production Systems

Zur Erlangung des akademischen Grades eines Doktor-Ingenieurs von der KIT-Fakultät für Informatik des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Julius Pfrommer

Tag der mündlichen Prüfung: 22. Juli 2019 Referent: Prof. Dr.-Ing. habil. Jürgen Beyerer Korreferent: Prof. Dr.-Ing. Michael Weyrich

**Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding the cover, pictures and graphs – is licensed under a Creative Commons Attribution-Share Alike 4.0 International License (CC BY-SA 4.0): https://creativecommons.org/licenses/by-sa/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2024 – Gedruckt auf FSC-zertifiziertem Papier

ISSN 1863-6489 ISBN 978-3-7315-1253-0 DOI 10.5445/KSP/1000152028

# **Kurzfassung**

Für automatisierte Produktionsanlagen gibt es einen fundamentalen Tradeoff zwischen Effizienz und Flexibilität.  Änderungen müssen aufwändig in einer Vielzahl von Systemen nachgezogen werden.

In dieser Dissertation wird ein Ansatz entwickelt, um eine automatische Anpassung des Verhaltens von Produktionsanlagen an wechselnde Aufträge und Rahmenbedingungen zu erreichen. Dabei kommt das Prinzip der Selbstorganisation durch verteilte Planung zum Einsatz. Die aufeinander aufbauenden Ergebnisse der Dissertation sind wie folgt:

Ausnutzung zusätzlicher Struktur in den Modellen skaliert der Ansatz auch auf große Szenarien.

Die Ergebnisse zusammenfassend werden Grundlagen für flexible automatisierte Produktionssysteme geschaffen – in einer Werkshalle, aber auch über Standorte und Organisationen verteilt – welche die ihnen innewohnenden Freiheitsgrade durch Planung zur Laufzeit und agentenbasierte Koordination gezielt einsetzen können. Der Bezug zur Praxis wird durch Anwendungsbeispiele hergestellt. Die Machbarkeit des Ansatzes wurde mit realen Maschinen im Rahmen des EU-Projekts SkillPro und in einer Simulationsumgebung mit weiteren Szenarien demonstriert.

# **Abstract**

There is a fundamental tradeoff between automation and flexibility in production systems. Large lot sizes can be produced efficiently with automated production systems. Many machines and equipment, like a 5-axis CNC mill, are in principle capable of producing many different kinds of parts. Similarly, (intra-) logistics systems exist for the automated transport and warehousing. Integrating these flexible components to an overall production system that is equally flexible has been prevented by the limits of automation technology in dealing with the ensuing complexity. Most production processes are rigid not only by way of the physical layout of machines and their integration, but also by the custom programming of the control logic for the integration of components to a production systems. Changes are time- and resource-expensive. This makes the production of small lot sizes of customized products economically challenging.

This thesis develops solutions for the automated adaptation of production systems based on self-organisation and distributed planning. The main results are the following:

duction processes. For this, the technique of Monte-Carlo Tree Search is extended. By exploiting algebraic structure in the models, the approach can be scaled to large scenarios.

In summary, the results of this work enable future production systems that are both efficient and adaptive. For this the components of the production system are enabled to use the degrees of freedom that are available to them. By the use of self-organization for the coordination, components of the overall system can react to changes in the system topology and external conditions. The approach was tested in application scenarios—in simulation and in physical production systems. For example as part of the EU-project SkillPro.

# **Acknowledgments**

Let me first thank Prof. Dr.-Ing. habil. Jürgen Beyerer for his mentorship, guidance and shared passion for agent-based and distributed systems throughout the development of this thesis. The time at the IES chair and the interactions were formative and invaluable. Furthermore I want to thank Prof. Dr.-Ing. Michael Weyrich for the discussions and his role as a reviewer of this dissertation.

Fraunhofer IOSB, and particularly the ILT department led by Thomas Usländer, was a great environment to conduct both scientific research leading up this thesis and to pursue high-impact engineering projects in the automation environment. The colleagues at IOSB and especially my former group leaders Miriam Schleipen and Ljiljana Stojanovic taught me the ropes of the craft of applied research. This was an experience from which I still profit immensely.

Besides the results on distributed planning for production control, a "byproduct" of this thesis was that it enabled me to contribute to the open62541 open source implementation of the OPC UA standard for industrial communication. This background ensures a clear technological path from the theoretical work to the application scenarios. The open62541 development team proved very inspirational and productive, so let me thank Florian Palm and Sten Grüner, Prof. Leon Urbas, Chris Iatrou, Stefan Profanter, and Andreas Ebner.

Another thank you goes to Joseph Warrington and Georg Schildbach for their guidance during my first steps into research during my time at the Automatic Control Laboratory at ETH Zurich.

Most importantly, this thesis would not have been possible without the encouragement and support of my family. *Ce thèse n'aurait pas été possible sans tout ton soutien, Ingrid!*

Karlsruhe, July 2019 *Julius Pfrommer*

# **Table of Contents**

# **List of Figures**

# **List of Tables**

# **List of Symbols**

#### **General Notation**

#### **General Sets**

### **Probability**

#### **Production System Model**

#### **Distributed Production System Model**

#### **Trace Theory**

### **Planning**

#### **Miscellaneous**

# **List of Acronyms**

# **1 Introduction**

*There has been a great deal of talk, much of it well founded, that the effect of science on economics and on the economy has not only been very large but that something like a second industrial revolution is impending. Illustrating this are the enormous advances in communications – physical and informational –, advances in automatization and in the domain of information and control, and finally, atomic energy.*

*John von Neumann [Neu55]*

## **1.1 Production and Logistics in a Global Economy**

After introducing the assembly line for the production of the Model T automobile in 1908, Henry Ford sought to make his company self-sufficient. For the supply of raw material, his Ford Motor Company bought 700,000 acres of forest, rubber plantations, iron and coal mines, and so on. The Ford River Rouge Complex near Detroit was designed to transform the incoming raw material into fully assembled cars. The manufacturing operations performed at River Rouge included coal coking, steel forging, sheet metal stamping, engine casting, lumber milling, tire making, the production of sheet glass from molten sand, and many more. All leading up to the final production step: the final car assembly in Ford's assembly line [Bri03]. Since Ford's only product at the time was the Model T, all production processes were highly specialized to maximize efficiency and reduce costs. This made the Model T the most affordable car at the time. Ford's competitor Chevrolet had a different approach. They used generic manufacturing equipment that could produce parts for several models at the same time. This also enabled frequent updates of the car models. Innovations of the 1920s, such as motors with electric starters that require no hand crank, let the Model T appear increasingly outdated. In 1927, Ford finally introduced a successor: the Model A. However, since the

#### 1 Introduction

production processes were tailored towards the Model T, the changeover proved difficult. The River Rouge site came to standstill for a duration of six months until production could be slowly resumed [Hou85]. In the following years, Ford gave up the model of the highly integrated manufacturing site. Today, like all automotive companies, Ford operates a range of production sites that are specialized on a range of parts (e.g. the internal combustion engine) that is used for several car models. And also the final assembly lines can switch between car models to adjust to changing customer demands. All production sites are connected with logistics networks and rely heavily on external material and component suppliers.

Ford's change from an integration manufacturing site to a network of interconnected production sites is exemplary for the evolution of many manufacturing industries. Initially, customized craft production is replaced by mass production, resulting in large efficiency gains and opening up new markets. When the markets are divided up, companies diversify their product portfolio to cater for individual customer groups. This leads to smaller order sizes, reducing the efficiency of the mass production approach. New methods have been developed to enable customization without loosing all the efficiency gains of mass production. Among the most popular ones are Lean Production [Ohn88] (also known as the Toyota Production System) and the Just in Sequence (JIS) inventory strategy [WS11]. They enable production sites to reduce the minimum order size that is still economical to produce. Sometimes lot sizes are reduced to the absolute minimum: a single customized product.

While automotive brands rely on a network of suppliers, these relationships are relatively stable. Building up the capacity to produce a specific part in high quantities in the expected quality takes time and investments in automated processes. This only makes economical sense when years of high demand are expected. For many retail goods, a similar division of labour happens in a complex supply-chain. But, these relationships can be established and dissolved literally overnight. Cheap long-range shipping, the reduction of tariffs and easy communication has enabled world-wide supply chains. We will take the example of global supply chains in the apparel industry [Ger99]. The supply-chain for a cotton shirt comprises the provider of raw cotton, spinning of the yarn, color dyeing of the yarn, weaving of the textile, cutting and sewing to make the shirt, design printing and stitching of brand logos. All of these production steps are typically executed by different companies in different countries. And the supply-chain is dynamically reconfigured for individual orders based on availability and price. A shirt bought today may have taken a wholly different way around the world than the same shirt bought a week earlier in the same store [Chr00]. Many apparel retailers even forego central warehousing for their stores. Instead, products are delivered directly from the last link in the supply chain to the store. So the retailers do not bind capital in stock for the entire season. And they can react within weeks to data showing good or bad sales of a specific product [CM15]. On the downside, most western apparel brands do not control their supply-chain and rely on sourcing agents from overseas.

The biggest sourcing agent is Li&Fung Limited. Operating out of Hong Kong, Li&Fung self-describes it's core business as "managing the supply chain for high volume, time sensitive goods" [FFW07]. Li&Fung owns no factory, warehouse or inventory and can still source nearly any retail good. Its database contains factories throughout Asia with their support for different manufacturing processes, available capacity and logistics options, as well as the availability and prices of raw material commodity components. Orchestrating the supply chain is a profitable business. In 2015, Li&Fung achieved a turnover of \$18.8 billion and a healthy \$2.2 billion profit [Lim15]. The leverage the sourcing agent has over the supply chain is seen increasingly critical by companies who rely on their services. In 2015, Wal-Mart announced a plan to reduce their reliance specifically on Li&Fung. Relying too much on a single provider had become a strategic weakness for the world's biggest retailer [WS15].

So why are dynamic reconfigurations of the supply-chain possible for the apparel industry, but not for automotive? In the apparel supply-chain, suppliers provide generic access to manufacturing processes that can be used for many different customers with little changeover costs and fast production ramp-up.

• Standardized commodity goods enable a high degree of automation. For example, the objective of textile plants is to run their power looms with as little downtime as possible. The automated equipment allows the configuration of different product types. For example, a Jacquard loom can configure different weaving patterns. The difference between product types can be entirely handled by the automated production system.

#### 1 Introduction

• On the other end of the spectrum, cutting and sewing of apparel is highly dependent on human labour. Lot-sizes are generally smaller and the type of product can change drastically between orders (e.g. switching from jeans to dress shirts). Handling of pliable textile material is difficult for automated equipment and requires custom machines and long changeover times. This makes automated equipment ineconomical as long as cheap human labour is available in overseas countries.

Globalization and the decentralization of supply chains lead to increased requirements for flexibility in production [Abe+06]. In practice, however, flexibility in production is a conflicting goal with efficiency. Automated systems provide the increased efficiency required for mass production. But they require considerable investment for the initial setup and the changeover between products. In recent years, many countries with a large industrial base have set up research programmes to renew industrial production with the increased use of information processing and communication technology. Among these programmes are "Industrie 4.0" in Germany [KWH13], "Made in China 2025" [Ken15] in China and "Industrie du Futur" [FD16] in France. A major part of these efforts is the creation of new automation technology that improves on the tradeoff between efficiency and flexibility in production.

## **1.2 The Structure of Automated Production Systems**

The vast majority of control systems for production systems is organized as a hierarchy. This is true both for discrete manufacturing and continuous production processes (e.g. chemicals, pharmaceuticals, beverages). Figure 1.1 shows the automation pyramid, a frame of reference for the hierarchical design of most automated systems in discrete manufacturing. Figure 1.2 depicts the typical automation hierarchy from the process industry (e.g. chemicals and pharmaceuticals). Decisions are made hierarchically and data is aggregated more and more as it is forwarded to the upper levels of the automation hierarchy. Hierarchical control follows the principle of subsidiarity, where upper levels make high-level decisions that are gradually refined as they are forwarded down the automation hierarchy. Subsidiarity is a necessary consequence of the fact that information

**Figure 1.1:** The automation hierarchy in discrete manufacturing. The IEC 62264 / ISA-95 standard does not include the enterprise level. It was added here to include interfaces to systems outside the shopfloor.

is aggregated when it moves up the automation hierarchy. The models used for decision making in the upper levels are more and more coarse and abstract away low-level details. But the low-level details have to be considered eventually. The lower levels in the automation hierarchy takes decisions from one level above and "fill the gaps".

In many systems, the control levels are tightly coupled. Changes to a component of a production system usually requires changes in many adjacent systems both vertically and horizontally. The subsystems are interwoven and implicit assumptions about adjacent systems are represented only in custom control code. This makes modifications to automated systems costly and time-intensive. Reducing this effort is the focus of an entire research community.

The remainder of this section gives an overview on the most important planning and optimization methods for decision-making on the different levels of the automation hierarchy. Whilst it is not possible to provide a complete enumeration, the examples from this chapter will set the frame for the modeling and planning techniques that are introduced later on.

**Figure 1.2:** Typical automation hierarchy of a chemical plant [Sko04].

#### **The Control Level**

Modern control theory and their implementation on computers can be traced back to work done at MIT in the 1940s. There, Norbert Wiener first coined the term "cybernetics" as the conjunction of control and communication [Wie48]. At the same time, project Whirlwind, conducted at Jay Forrester's Servomechanisms Laboratory, developed digital "feedback control" for numerically controlled (NC) manufacturing processes [Rei91]. Both modern Programmable Logic Controllers (PLCs) and the application of feedback control methods on digital computers are descendant from this work in a direct lineage. The fundamental difference between the two lies in the type of decisions they have to make. Programmable Logic Controllers (PLC) generally are driven by a state machine with discrete transitions or events. Feedback (optimal) control is mostly concerned with physical systems with continuous dynamics.

**Figure 1.3:** Example for a PLC program in ladder logic

Programmable Logic Controllers Programmable Logic Controllers (PLC, see [Joh87; Wal12]) are directly interfaced with a physical process via sensors and actuators. In the automation of discrete manufacturing, PLC couple the physical system with digital control and communication. This coupling requires custom program code to accomodate for the specific details of the physical system and its intended functionality. As the control level often deals with safety-critical functionality, hard bounds on realtime reactiveness have to be guaranteed. The IEC-61131 languages [Com93] are standard for PLC programming and mandate a programming style that is idiomatic to industrial controllers. Figure 1.3 shows an example for one of the IEC-61131 languages, ladder logic, which is directly descendent from the analog circuits that were originally used for industrial control before the PLC.

The logic coded into a PLC is mostly reactive. Sensor inputs are read and translated into actuation commands. This loop is repeated at a fast pace (many hundred Hertz) for realtime operations. Lengthy planning procedures do not usually fit into the constraints of the control loop in a PLC. The IEC-61499 standard is intended as a modernization of IEC-61131 [Vya11]. It adds an event-based control flow to the strictly cyclic operations of previous PLC generations. But even with IC-61499, PLC-based control is mostly reactive. Computationally expensive planning is generally not performed in a safetycritical control environment.

Automated production systems typically are integrated from components that come with their own control hardware. The integration requires communication between individual controllers and custom control software to react to cyclically transmitted status messages and events. On the control level, traditional fieldbuses are still common today [Zur14]. Fieldbuses are even mandatory to use if safety-critical functionality relies on digital communication between controllers or between a PLCs and field devices with sensors and actuators. The programming of PLCs is often finished on-site as part of the integration of system components. As many automation systems are custom solutions, code reuse for PLCs is difficult even if systems are built from standard components. This leads to a tight coupling that also increases the time required to make changes in an existing system.

Feedback Control The canonical definition of an optimal control problem is as follows [Lib11]:

$$
\dot{x} = f(t, x, u), \quad x(t\_0) = x\_0 \tag{1.1}
$$

The system dynamics is described by an ordinary differential equation (ODE) *f*. The system state at time *t* is x(*t*) ∈ R *<sup>n</sup>* and its evolution depends on the control input u(*t*) ∈ R *<sup>m</sup>*. The initial condition is given by x0. A cost functional for the state evolution assigns costs to the state and control effort between times *t*<sup>0</sup> and *t<sup>f</sup>* and an additional terminal cost on the final state x(*t<sup>f</sup>* ).

$$C(\mathfrak{u}) = \int\_{t\_0}^{t\_f} L(t, \mathfrak{x}(t), \mathfrak{u}(t))dt + K(t\_f, x\_f) \tag{1.2}$$

The problem of optimal control is to compute u <sup>∗</sup> = arg min<sup>u</sup> *C*(u). In this general framework, computing u ∗ is a variational problem as u is a (vectorial) function over time [Lue69]. A popular approach is to discretize the time domain *T* = {*t*0*, t*1*, . . . , t<sup>f</sup>* } so that u <sup>∗</sup> problem of finding the sequence of u that minimizes *C*. Applying the resulting u blindly until *t<sup>f</sup>* is called *open-loop* control. Repeating the optimization after every time period with updated state information is called Model Predictive Control (MPC) or Receding Horizon Control [ML99; Mac02].

Traditionally, feedback control has be implemented as analog electrical circuits. But these are restricted to relatively simplistic solutions, such as PID controllers [Ben93]. With the increase in computational power available in control devices, Model Predictive Control (MPC) has become possible for many application. Here, the control problem is stated as an optimization problem that is solved repetitively at a high frequency to incorporate sensor measurements for "feedback" control. Optimal control as an optimization problem originates from the association of dynamical system with control input with

## **The Operations Level**

or a buildup of residual oil from the metal coating in the stamp tool. It is often up to a skilled process expert to adjust the process parameters at runtime to ensure the required performance.

### **The Plant Management Level**

On the plant management level, an entire shopfloor is considered. Typically the planning horizon is between several hours and several days of operation [She03].

### **The Enterprise Level**

Enterprise Resource Planning Enterprise Resource Planning (ERP) describes a class of software systems to assist enterprise-wide management [Jac+07]. The

ERP products with the highest market share are SAP ERP (previously SAP/R3) and the Oracle E-Business Suite [Gar18]. In many aspects, ERP functionality mirrors tasks performed at the plant management level. But the timeframes are generally much longer and several production and warehousing sites are jointly considered. Enterprise-wide optimization (EWO) aims at optimizing the operations of supply, manufacturing and distribution activities of a company to reduce overall costs and inventories [Gro05].

Supply-Chain Management Supply-Chain Management (SCM) [Ali05; GF08] is concerned with production scenarios that include several layers of suppliers. SCM is most common in industries where suppliers are not delivering specialized parts instead of commodity products. The integration with suppliers is often very tight. This enables the reduction of buffer storage at the production site by the use of just-in-time and just-in-sequence delivery. The so-called bullwhip effect [LPW97] describes how small fluctuations in customer demand lead to large fluctuations in demands at suppliers that are several tiers removed. A big motivator for digitalisation and information sharing in the supply-chain is the reduction of the bullwhip effect by enabling better forecasts for the suppliers.

## **1.3 Approaches for Flexible Production Systems**

This section discusses approaches to render automated production system flexible. The possibility of flexible production is one of the driving motivations for Industrie 4.0 [Wey+14]. This thesis has a scope on automated production. Organizational methods that focus on the human element in production, such as Lean Manufacturing [Ohn88], are therefore not considered in depth. Several authors have developed frameworks to characterize flexibility in production systems [Bro+84; BS88; GG89; SS90; Ger93; DT98; Wie+07]

The scientific literature uses specific terms to describe flexibility in a production context. See Figure 1.4 for a common nomenclature by Wiendahl [Wie+07]. Even though specific terms exist, the term flexibility is deliberately used with its colloquial meaning in this thesis: The models and planning algorithms developed in this thesis apply to all level in the automation hierarchy. The specific terms for flexibility from the scientific literature are mostly tied to one level of the

**Figure 1.4:** Classes of manufacturing changeability from [Wie+07].

automation hierarchy. We aim to avoid misunderstandings by the specific terms outside of their commonly understood definition.

**Service-Oriented Production Systems** The principle of service-orientation is used in computer science to develop system architectures where components are loosely coupled [Mac+06]. A specific service provider can be exchanged as long as the interfaces for interaction remain identical and the underlying functionality is still provided. Discovery mechanisms are used to find and select appropriate service providers. In the context of Industrie 4.0, service-orientation is regarded as an enabler for future control system architectures that dissolve the classical automation hierarchy. See Figure 1.5 for a popular depiction. The DIN SPEC 16593-1 standard [DIN18] defines a reference model with basic principles for service-based architectures in the context of Industrie 4.0. This is the common basis for technical realisations to the vision from Figure 1.5.

**Figure 1.5:** Decomposition of the automation hierarchy with distributed services [Mes13; Mon14].

A range of research projects has translated service-orientation to production control. The authors from the SOCRADES project [JS05; De +08; Cân+11] and Shen et al. [She+07] develop a service-oriented manufacturing system architecture where semantic technologies are used to match possible providers of functionality in a manufacturing system. Loskyll et al. [Los+11; Los+12] expand the concept of semantic service discovery to the parameterization and orchestration of services. For this, they develop a domain-specific ontology for the use in semantic reasoning tools. Puttonen et al. [PLM13] describe a set of specialized web services for composing and invoking semantically enriched automated procedures in a manufacturing setting. They also present an algorithm to identify the steps required to reach a predefined goal state. Dürkop et al. [Dür+14] discuss the use of service-oriented architectures in reconfigurable manufacturing systems (see Figure 1.4 and the technical challenges that need to be overcome. [SZW17] use model-based approach for the service development and a modular architecture to reduce the complexity of service-based production systems. [LV15] combine service-oriented manufacturing control with a multi-agent architecture. The Smart Factory Web testbed in the Industrial Internet Consortium (IIC) uses web services for planning and control in global supply chains [Jun+17].

**Agent-Based Production Systems** An even more radical departure from the classical automation hierarchy is investigated with agent-based distributed control of production systems. Agent-based systems in manufacturing and logistics are the topic of a dedicated research community that has been active since the 1980s [DP87; LS92]. The survey papers [MVK06; LK08; Lei09; LMV13; LK15] give an account on the history of agent-based control and an overview on the focus of current work. Notably, the IEEE-IES Technical Committee on Industrial Agents (TC-IA¹ brings together researchers on an international level.

Software frameworks have been developed to assist the development of agentbased systems. For example the well-known JADE project [BCG07]. The frameworks for software agents, however, do not provide abstractions specifically for the production domain and are used for the development of distributed software systems in general.

The core challenge of agent-based control is the coordination of individual decision making across agents. The remainder of this paragraph discusses the most common approaches. A common coordination mechanisms for agent-based control is negotiation [ZR89]. The Contract Net Protocol (CNP) [Smi80] replaces the market with a negotiation scheme in order to decompose and distribute tasks between agents. The CNP has been applied for agent-based manufacturing systems in a range of research projects and industrial installations [Par87; LL94; SKB97; Oue+99]. While the CNP is mostly used for greedy decision making, other authors have integrated scheduling theory with agent-based control [SWH06; Agn+14; Bad11] Other coordination mechanisms are nature-inspired and derived from the behavior of animals [XL08].

Holonic production control is a special case of agent-based control. The term *holon*, originally coined in [Koe68], refers to systems made up from components that encapsulates both physical assets and virtual functionality [GLK98; Fis99; MB00]. The key idea is that the system components are themselves holons. This goes beyond the usual system-of-systems approach, as holons are self-similar in the sense that the structure and functionality of the constituent parts is governed by the same principles as their parent. Taken to its extreme, this self-similarity in manufacturing systems has led to the concept of the fractal factory [War93].

¹https://tcia.ieee-ies.org/

**Figure 1.6:** Building blocks of a holonic manufacturing system according to the PROSA reference architecture [Van+98].

The PROSA project has proposed a architecture reference architecture for holonic manufacturing systems [Van+98]. See Figure 1.6 for the building blocks defined by PROSA.

What is currently lacking in the field of agent-based and holonic manufacturing control are widely used benchmark scenarios to quantitatively compare the proposed coordination mechanisms. Compared to other scientific fields, this has led to many competing approaches without a clear winner and uncertainty on how well the different approaches can cope with aspects outside of their original scope. For example if unforeseen events are introduced in a stochastic environment. For practitioners, this has led to an overwhelming range of choices. For researchers, years of effort have so far not amalgamated into a unified theory of agent-based production control.

### **Plug and Produce**

The idea of Plug & Produce is derived from plug-and-play functionality known from the USB interface for computer hardware. There, well-known device classes with standardized functionality remove the need for custom software drivers for the hardware integration. Arai et al. [Ara+00] first translated plug-an-play to the production domain and coined the term Plug & Produce. Onori et al. [Ono+12] use the concept for a self-configuring assembly system at the shop-floor level.

The integration of machines and equipment with Plug & Produce encompasses the following aspects: First of all, basic connectivity is established for an existing (industrial) communication infrastructure [Dür+12; Rei+10]. Second, the new component announces its presence to a central controller or directly to the adjacent components with a discovery mechanism [Pro+17]. Third, in production, there exists a wide range of machines and equipment. This heterogeneity cannot be reduced to a small number of devices classes. A way to enable Plug & Work scenarios in the face of device heterogeneity is the use of self-descriptions languages for the integration [OHN14; Sch+15a]. The fourth and most challenging aspect is the functional integration. Lepuschitz et al. [Lep+11] show reconfiguration of manufacturing resources based on distributed IEC 61499 function blocks and a semantic description of the manufacturing setting and the expected behavior.

Many of the published Plug & Produce implementations use a dedicated interconnector module that acts as a facade for manufacturing equipment and provides a uniform interface and that generates low-level commands for the underlying device [NWS07; Dor+17].

One approach for the functional integration in Plug & Produce is the modeling of the skills of technical equipment. For this, see the review of the state of the art in Chapter 5.

## **1.4 The Missing Hierarchy of Production Theories**

Scientists and engineers use models on a level of abstraction that is the most useful for the phenomena under investigation [Gie04]. It often occurs that a more accurate model is available in principle. But working with an increased level of accuracy would overburden the analysis with unnecessary complexity. For example, an electrical engineer laying out the power grid of a city will not use Maxwell's Equation for a power-flow study. Many technical fields have arranged these model approaches (theories) in a hierarchy. This hierarchy has evolved over time and – in the natural sciences – its development is closely related to the process of scientific discovery [Kuh62; Car84]. As an example, Figure 1.7 shows how the model hierarchy established in the field of optics . If some phenomena cannot be

**Figure 1.7:** The model hierarchy in optics. (Reproduced from http://www. argmin.net/2018/01/25/optics.)

explained one can resort to a more detailed (and computationally or analytically more expensive) model until the first principles from Maxwell's Equations and quantum physics are reached [MW59].

Different academic fields have produced "Theories of Production". For example economics [Sch34; Dan66; She71], business administration [Sch86; Fär88] and production management [Dyc06]. Around the year 2000, prominent authors have called for a unified theory of production that provides a common foundation that integrates existing results [Dyc03; Sch04; WNH10]. Recent years have seen a range of proposals to fulfill this need [NW10; Sch+11; Sch+15b]. The authors of [Sch+17] provide a comprehensive review and classification of such production theories. But the cited work remains mostly conceptual and does not model the detailed control of automated systems on the lower levels of the control hierarchy.

By contrasting the model hierarchies in other academic fields with the control hierarchy in automation, one could aim for a high-fidelity model at the bottom layers that is abstracted more and more as we go up in the control hierarchy. While that is the case, we want to stress that the underlying modeling principles should stay the same for all levels of the control hierarchy.

The work on a theory of production is not only of academic interest and also relevant to practitioners that have to cope with the increasing complexity of production systems. Currently, each layer of the control hierarchy works with dedicated system models that are tailored for the tasks at hand. But the models are so task-specific that they become mutually incompatible. This is problematic at vertical as well as horizontal interfaces of the control hierarchy. Assumptions about the behavior of adjacent components are implicit in the custom rules for the interaction between subsystems and encoded in custom control program code. The lack of a common core to translate between subsystems leads to constant manual effort for system integrators and the loss of flexibility due to the task-specific hard coupling of components.

## **1.5 Scientific Contributions and Thesis Organization**

The original title at the very beginning of the work leading to this thesis was "Agent-based production control using the paradigm of self-organisation in the context of Industrie 4.0". The path to the final thesis led through scientific fields that are not directly associated with that original title. For example Probabilistic Graphical Models, Convex Optimisation, Description Logics, Reinforcement Learning, and many more. Some of these detours did not pay off as intended. But many did. While not always visible, the thesis has stayed true to the original topic and the results and techniques from different fields finally tie to together into one body of work.

The thesis is organized into six chapters. See Figure 1.8 for an overview and recommended reading order. The scientific contributions are closely related to the thesis structure. In general, prior results are summarized in dedicated "background" sections. Beyond the background sections, all results with an explicit reference to the literature were developed as part of the thesis development. It follows a summary of the key developments.

**Chapter 2: A Model of Concurrent Production Systems** A novel model representation for production systems is introduced. The goal of the model is to provide a uniform representation of both discrete and continuous production

**Figure 1.8:** Structure and recommended reading order of the thesis.

processes on all levels of the automation hierarchy. To achieve this, the model combines the following aspects:

The following four chapters aim to show that the model representation not only presents a common core for higher-level models but is practical to use in flexible production systems. For this, a tailored planning algorithm is developed,

further scaling of planning is achieved by exploiting additional model structure, the planning algorithm is extended to distributed planning by independent agents, and finally models of the skills of machines and equipment are used as highlevel abstraction from which executable low-level representations are derived for runtime control.

#### **Chapter 3: Simulation-Based Planning for Concurrent Production Systems**

The production system model implies a planning problem: maximizing the expected reward that is generated. Existing planning algorithms from the production domain do not apply to the model from Chapter 2. They make limiting assumptions on the planning problem structure that no longer apply. To solve the planning problem, an algorithm based on Monte-Carlo Tree Search (MCTS) is developed.

**Chapter 4: Distributed Planning for Self-Organizing Production Systems** The planning complexity depends on the number of individual components that partake in the system. By reducing the planning scope to a subset of the system components, it is much easier to explore the solution space and converge to good solutions. We develop an algorithm that combines Monte-Carlo Tree Search with Message Passing approaches known from Belief Propagation. Thereby, the system can be compartmentalized into individual agents who are coordinating their actions. The agent coordination mechanism is shown to improve the planning solution quality compared to uncoordinated individual planning. The same implementation of our novel planning algorithm can be used the solve hitherto separate planning and optimization concerns from all levels of the automation hierarchy.

**Chapter 5: Modeling of Production Skills** In order to achieve flexibility in automation with runtime planning, an accurate system model is required. It can be quite resource-intensive to keep the physical system and its model representation synchronized. Especially if new manufacturing operations are frequently introduced in a flexible production environment. To reduce this effort, higher-level descriptions are developed from which low-level representations for runtime control are generated. We introduce a formal model for the technical skills of components in an automated system based on semantic modeling and deductive inference based on second-order logic.

The modeled skills are used to generate executable action-representations for specific operations. The generation of executable actions take in as input the capabilities of the system components, the topological layout of the system and a description of the requested operation. The capability model is therefore used as a high-level descriptive language that can be compiled to executable actions for the participating system components. The generated action representations contain all preconditions and effects required for detailed planning.

Already published results that were created in preparation for this thesis are:

The following publications are our technical reports or publications in adjacent fields. They contributed to the thesis indirectly by using similar techniques or by working on implementation technologies in industrial automation (e.g. using technologies such as OPC UA and AutomationML) that are relevant for bringing the results of this thesis into practice.

# **2 A Model of Concurrent Production Systems**

*The value of [Lagrange's book "Mécanique Analytique"] consists in the exposition of a general method by which every mechanical question may be stated in a single algebraic equation. The entire history of any mechanical system, as for example, the solar system, may thus be condensed into a single sentence.*

*Robert S. Woodward [Woo95]*

The chapter introduces a model for production systems that combines continuous and discrete system dynamics with concurrency, i.e. parallel operations and the synchronization of system components. The model is intended as the basis for a "theory of production systems". For this, the model needs to able to represent the system dynamics on all levels of the control hierarchy. Examples from different levels of the control hierarchy are used for demonstration and to substantiate this claim. The core postulate of this chapter is the following:

> *The same set of modeling principles can represent the relevant properties of production systems on all levels of the control hierarchy.*

## **2.1 State, Actions and Action Sequences**

**Definition 2.1.** *A system is a set of components c* ∈ *C.*

Examples for components on a manufacturing shopfloor are machines and logistics equipment, such as forklifts. Depending on the level of abstraction of the model, components can also represent parts inside a machine, as well as entire production plants and warehouses. Assume for now that components do not contain other components in a hierarchy.

**Definition 2.2.** *Each component c has a state s* ∈ *Sc. The set S<sup>c</sup> contains all possible states of the component c.*

The global state space *S* = ×*c*∈*<sup>C</sup> S<sup>c</sup>* is the cartesian product of the possible states for every component. A global state is a vector s ∈ *S* with elements *sc*. The state space of a subset of the components *Q* ⊆ *C* is *S<sup>Q</sup>* = ×*c*∈*QSc*. The projection Π*Q*(s) = (*s<sup>c</sup>* : *c* ∈ *Q*) extracts the state of the components *Q* from the global state. In general, a subscript denoting a set of components indicates projection and s*<sup>Q</sup>* = Π*Q*(s). The inverse projection is Π −1 *<sup>Q</sup>* (s*Q*) = {u ∈ *S* : Π*Q*(u) = s*<sup>Q</sup>* }.

**Example 2.1.** Consider the ABB IRB140 robotic manipulator from Figure 2.1. The green sphere represents the Tool Center Point (TCP). The TCP has six degrees of freedom (three each for translation and rotation). But some positions are not reachable due to the physical constraints of the manipulator. The constraints are encoded in the set of reachable positions Ψ ⊂ R 6 . In addition, the IRB140 can be fitted with different tools. In this example, the possible tools are Υ = {gripper*,* welder*,* drill*,* none}. The overall state space of the manipulator is *S*irb140 = Ψ × Υ.

**Figure 2.1:** IRB140 robotic manipulator with inverse kinematics control (Coppelia V-Rep).

Many manufacturing operations define trajectories where the manipulator starts and finishes at fixed locations. As an alternative to modeling a continuous state space, one can limit the possible positions to a discrete set of predefined positions. A continuous configuration space for the position is however better suited to model the physical movement dynamics of the robot.

Components can change their state over time. The following nomenclature is taken from [Fuj98]: *Physical time* refers to time in the physical system that is represented by the model. *Simulation time* refers to a time representation, for example a real value that corresponds to a physical clock by scaling and an epoch-date for the origin. *Wallclock time* refers to the time when the simulation is executed. In distributed systems, it is generally impossible to assign an absolute order to events [Lam78]. We make the simplifying assumption that clocks are synchronized to absolute precision. So the physical time of the system components is always identical. In the model representation, however, the simulation time can differ between components. That is, components can evolve their state independently from one another until synchronization forces their simulation time to coincide again. If not stated otherwise, time refers to simulation time in the remainder of the text.

**Definition 2.3.** *The time-indexed state of a component c is*

$$
\sigma = (s, t) \in \Sigma\_c, \quad \Sigma\_c = S\_c \times \mathbb{R} \,. \tag{2.1}
$$

*It indicates the state of the component s at time t. The time t represents the offset from some epoch-date in seconds.*

The time-indexed state of the entire system is σ ∈ Σ for Σ = ×*c*∈*<sup>C</sup>* Σ*c*. In the context of a system state σ, the time-indexed state of a component *c* is referred to as *σ<sup>c</sup>* = (*sc, tc*). Again, subscript-based notation for projection refers to the definitions from the surrounding context: Σ*<sup>Q</sup>* denotes the joint time-indexed state-space for a subset of the components *Q* ⊆ *C* and σ*<sup>Q</sup>* = Π*Q*(σ) is the projection of a global time-indexed state into Σ*Q*.

Components change their time-indexed state by executing actions. In the model, actions skip the time-indexed state ahead to the time after the execution. Components have no well-defined state during the execution of an action. If the model is linked to a physical instance of the system, then the relation between simulation time and physical time is as follows. If the time *t<sup>c</sup>* of some component *c* is earlier or equal to the physical time, then the component has been idle since *t<sup>c</sup>* and is immediately available. If *t<sup>c</sup>* is later than the physical time, then the component is occupied with the execution of one or more actions until *tc*.

**Definition 2.4.** *Actions are possible state transitions, represented as a tuple*

$$a = \left(C\_a, \bar{\Sigma}\_a, \mathbf{e}\_a, d\_a\right). \tag{2.2}$$

*The tuple describes the action via its participating components, preconditions and effects. It consists of*

The set of feasible time-indexed initial states Σ¯ *<sup>a</sup>* encodes the preconditions of the action *a* for the participating components. The shorthand notation for the state of the participating components is σ*<sup>a</sup>* = Π*C<sup>a</sup>* (σ). The feasible global initial states for the action *a* are Σ*<sup>a</sup>* = {σ ∈ Σ : σ*<sup>a</sup>* ∈ Σ¯ *<sup>a</sup>*}. Every action is an operator on the time-indexed global state *a* : Σ*<sup>a</sup>* → Σ.

From the initial global state σ it follows that the earliest possible starting time is *t* start *a* (σ*a*) = max*c*∈*C<sup>a</sup> t<sup>c</sup>* for the action *a*. Let σ <sup>0</sup> = *a*(σ) denote the global state following the execution. The new global time-indexed state is then comprised of elements

$$\sigma\_c' = \begin{cases} \left( e\_{a,c}(\sigma\_a), \ t\_a^{\text{start}}(\sigma\_a) + d\_a(\sigma\_a) \right), & \text{if } c \in C\_a \\ \sigma\_c, & \text{else} \end{cases} \tag{2.3}$$

Equation 2.3 defines the operational semantics of all actions. It implies the *Markov assumption*: The outcome of an action depends only on the previous system state. Components that are not participating in an action do not partake in the preconditions and their state is left unchanged as well.

Equation 2.3 defines how components with different simulation times are synchronized: The participating component with the highest simulation time defines the starting time of the action. The other participating components idle until they join the execution.¹ The feasible initial states Σ*<sup>a</sup>* can also encode preconditions for the component timing. Consider a physical system whose dynamics is described by a differential equation. Such systems typically cannot idle without any effect on their state. To adequately cover components that cannot idle in the sense of Equation 2.3, the preconditions of Σ*<sup>a</sup>* can require identical simulation times for all participating components ∀σ*<sup>a</sup>* ∈ Σ*a,* ∀(*c, c*<sup>0</sup> ) ∈ *C<sup>a</sup>* × *Ca, t<sup>c</sup>* = *t<sup>c</sup>* <sup>0</sup> . Alternatively, as e*<sup>a</sup>* takes the time-indexed state of all participating components as input, the state evolution of the individual components until the starting time of *a* can be modeled as part of the effect function.

In discrete manufacturing, the major concern is the movement and transformation of products within the system. Many modeling approaches represent products as objects with individual lifecycles [Sal+10] or even as independent agents who negotiate and make independent decisions [KBT17]. To model products as individual objects, they could be represented as system components *c*. But, in this text, we follow a different approach that enables algorithmic improvements for planning later on.

Let *ρ* ∈ *P* denote the set of different product types. A product type not only represents marketable output, but also raw material and semi-finished work-pieces that occur between production steps. Products are not considered individually. So products of the same product type are indistinguishable. Instead of representing each component individually, it suffices to track the number of products of each product type present at every component. A product is always contained in some component. Take for example a workpiece mounted inside a machine, a crate of material sitting on a forklift, or a finished product stored in a high bay warehouse.

¹This corresponds to the use of the popular Max-Plus algebra to describe the time-evolution of discrete event systems [Bac+92].

**Definition 2.5.** *Some components can physically contain products. In that case, the component state decomposes into the component configuration ξ* ∈ Ξ*<sup>c</sup> and a vector* p *for the number of contained products for each product type.*

$$s = (\xi, \mathfrak{p}) \in S\_c, \quad S\_c \subseteq \left(\Xi\_c \times \mathbb{N}\_0^{|P|}\right) \dots$$

For convenience, we denote a vector with just a single product of type *ρ* as the basis-vector ν*ρ*. The number of products of type *ρ* in component *c* is (*pc*)*ρ*. A component *c* with no contained products has p*<sup>c</sup>* = 0, where 0 denotes the null-vector of appropriate dimensionality.

**Example 2.2.** This example introduces a minimal manufacturing scenario that will be used throughout this text. The scenario, shown in Figure 2.2, is comprised of a machine tool (mt) that mills piston rods out of steel bars, a lattice box (box) and a robotic manipulator that packages the piston rods for transport (manip).

**Figure 2.2:** Minimal scenario for discrete manufacturing

Four actions have been defined for the scenario. The actions produce and package have only one participating components, the machine tool and the manipulator respectively. The actions put and take require the participation of two components to model the transition of a product between them. We do not consider single piston rods, but only orders of 100 parts with the product type order. Every component can only hold a single order-product at once. Therefore, once the machine tool has completed an order, the parts need to be put into the lattice box before the next order can be produced. The manipulator then has to take the parts of the previous order out of the lattice box before the machine tool can put in the next. Products "appear" and "disappear" when they leave the scope of the modeled system.

The complete definition of the action put is as follows:

$$\mathbf{\bullet} \cdot C\_{\mathbf{put}} = \{\mathfrak{mt}, \mathbf{box}\}$$

$$\begin{aligned} \bullet \; \bar{\Sigma}\_{\mathtt{put}} = \{ \sigma\_{\mathtt{put}} \in \Sigma\_{C\_{\mathtt{put}}} : p\_{\mathtt{nt}} = \nu\_{\mathtt{order}}, \; p\_{\mathtt{box}} = 0 \} \end{aligned} $$

$$\mathbf{\bullet} \,\, e\_{\mathrm{put}}(\sigma\_{\mathrm{put}}) = ( (\xi\_{\mathrm{mt}}, \mathbf{0}), (\xi\_{\mathrm{box}}, \nu\_{\mathrm{order}}) )$$

$$\mathbf{\bullet} \cdot d\_{\text{put}}(\sigma\_{\text{put}}) = (5s, 5s)^2$$

Actions can be chained to form action sequences. The set of actions *A* is the alphabet for the free monoid (Kleene Star, [HU79]) *A*<sup>∗</sup> which contains all words (sequences) of finite length. The empty sequence is written as *ε*. The sequence elements *w k* are indexed as w = *w* <sup>1</sup>*w* 2 *. . . w*|w<sup>|</sup> .

**Definition 2.6.** *An action sequence* w ∈ *A*<sup>∗</sup> *is itself an action resulting from the composition of the constituent actions.*

For a valid initial state σ ∈ Σ<sup>w</sup> (defined in the next paragraph), the resulting state after the execution of w is w(σ) = (*w* <sup>|</sup>w<sup>|</sup> ◦ · · · ◦ *w* 1 )(σ) according to the dynamics of action execution from Equation (2.3). The multiplication notation for sequences is preferred to the ◦-notation for function composition to keep the notation light, to write the sequence elements in-order, and because we will take on an algebraic perspective on composition later on. Note that, due to the way concurrency is represented with the time-indexing of component states, an action at a later sequence index may actually start at an earlier simulation time than one of its predecessors in the sequence. But this is possible only if the two actions do not share participating components.

The subsequence of the first *k* elements is w:*<sup>k</sup>* = Q*<sup>k</sup> <sup>l</sup>*=1 *w l* . The subsequence starting at the *k*th element is w*k*: = Q<sup>|</sup>w<sup>|</sup> *<sup>l</sup>*=*<sup>k</sup> w l* . The domain of w as an operator is Σ<sup>w</sup> = σ ∈ Σ : ∀*k* ∈ {1*, . . . ,* |w|}*,* w:*k*−<sup>1</sup> (σ) ∈ Σ*w<sup>k</sup>* . This ensures that the preconditions of all actions are satisfied when the sequence is executed in order.

The composition of action operators to a sequence is always possible. But the resulting sequence might be infeasible with an empty domain Σ<sup>w</sup> = ∅.

The set of all action sequences, from now on denoted as *W* = *A*<sup>∗</sup> , implies a tree-graph. In the context of an initial system state σ, infeasible sequences are removed from the tree. The feasible sequences starting at σ are *W*<sup>σ</sup> = {w ∈ *A*<sup>∗</sup> : σ ∈ Σw}.  After the first action *a* has been executed, the system states becomes *a*(σ) and the sequence tree becomes *a*(*W*) = {w : *a*w ∈ *W*}.

**Example 2.3.** This extends the scenario from Example 2.2. Assume an initial system state σ where no component contains any products. The four defined actions with their preconditions and effects yield a sequence tree *W*<sup>σ</sup> of feasible sequences. All feasible sequences up to five actions are shown in Figure 2.3.

**Figure 2.3:** Sequence tree for Example 2.2.

Note that the order of the actions in the sequence do not require that the start times of the actions have the same ordering. Consider the action sequence (produce put produce take). The last two actions produce and take can start at the same time. It is even possible that an action occurs later in a sequence but starts before a preceding action in terms of absolute time.

**Table 2.1:** Evolution of the system state for the sequence produce*,* put*,* produce*,* take. The component state *s* is simplified and describes only the number of products currently at the component. All actions are assumed to have a duration of five time units.

Action sequences can be the result of a planning procedure that are executed only once for the current situation. But sequences can also be used as reusable macro-actions. Many manufacturing systems are organized in lines, within which work-pieces undergo a fixed sequence of actions. In Example 2.3, the sequence (produce put take package) could be such a macro action. System parts with little flexibility can thus be modeled with comparatively few macro-actions. These can be integrated seamlessly with more fine-grained actions where more behavioral flexibility is required.

## **2.2 Parameterized Actions**

The action definition from Section 2.1 has been accompanied by examples from discrete manufacturing. In the process industry (e.g. chemicals, pharmaceuticals, food and beverages), many decisions have to be made on a continuous domain. In principle, the set of actions *A* could represent continuous decisions with an uncountably infinite set of actions. But then we could no longer maintain explicit representations of every action in computer-based simulations. Instead, we allow actions to be parameterized.

**Definition 2.7.** *The preconditions, effects and durations of a parameterized action a depend on the choice of action parameter θ* ∈ Θ*a.*

$$a^{\theta} = (C\_a, \bar{\Sigma}\_a^{\theta}, \mathbf{e}\_a^{\theta}, d\_a^{\theta}) \tag{2.4}$$

There are no particular restrictions on possible parameter spaces Θ*a*. Of course, parameters can also be vectorial or sets, even though we use scalar notation for parameters in general. Actions *a* that do not define parameters have Θ*<sup>a</sup>* = ∅. In that case, the parameter can be omitted in the notation. The following list gives examples for parameters on different scales of measure and sizes of the parameter space.

Nominal and Ordinal Parameters In the simplest case, parameterized actions simply group actions that, in some sense, belong together. For example the action paint with categorical parameters Θpaint = {blue*,* green}. There can also be a natural order among the parameters on an ordinal scale, such as {cold*,* warm*,* hot}.

**Example 2.4.** Take the example of a storage tank in a process control setting. The tank can be filled with liquid and drained afterward. The fill level is controlled by a pump and a valve at the bottom of the tank.

**Figure 2.4:** Piping and instrumentation diagram (P&ID) of a storage tank. (Reproduced with permission from https://commons.wikimedia. org/wiki/File:Pump\_with\_tank\_pid\_en.svg.)

Consider the action drain acting on the tank. Let *v* ∈ [0*,* 1] be the control value for the valve and *τ* the duration of the action with *θ*drain = (*v, τ* ). The valve is closed for *v* = 0 and fully open for *v* = 1. In addition, the liquid flow depends on the pressure at the valve and hence the fill level described by the state of the tank *s*tank ∈ R+. The flow of the (incompressible) fluid out of the tank can then be described by a differential equation *s*˙tank = *g*(*s*tank*, v*) according to Bernoulli's principle [Ber38] and *s*tank(*t*0) is known from the initial state σdrain. The effect of the action drain on the tank is

$$e^{(v,\tau)}\_{\mathbf{d}\mathbf{r}\mathbf{a}\mathbf{in}}(\sigma\_{\mathbf{d}\mathbf{r}\mathbf{a}\mathbf{in}})\_{\mathbf{t}\mathbf{a}\mathbf{k}} = s\_{\mathbf{t}\mathbf{a}\mathbf{k}}(t\_0) + \int\_{t=t\_0}^{t\_0+\tau} g(s\_{\mathbf{t}\mathbf{a}\mathbf{k}}(t), v)dt\ .$$

With the introduction of parametric actions, problems from control theory can be represented in the model. In Model Predictive Control (MPC) [ML99], a dynamical system is approximated by discretizing the time domain. The problem of optimal control is then posed as an explicit optimization problem that selects control values for every time period. By taking the (continuous) control decisions as action parameters, optimal control problems can be represented with a single action that encompasses all system components. Control with mixed discretecontinuous control values can be represented by either a more complex parameterspace Θ or by representing discrete choices with different actions. But only with several actions can aspects of concurrency be represented, where the system components are synchronized via joint participation in an action and evolve their state independently otherwise.

So far, the model presented in this chapter unifies the treatment of discrete events, concurrency and continuous system dynamics.

## **2.3 Uncertainty and Observations**

Until now, it was implicitly assumed that actions are deterministic. But virtually no production system actually is. There is always an interaction with a stochastic environment, from logistics influenced by weather conditions to the human operator returning late from a break. In addition, even fully automated manufacturing processes are inherently stochastic. This is illustrated by the fact that very few production processes achieve zero-defects production.² It is desirable to represent this uncertainty also in the model used for control or planning. With an explicit representation of the uncertainty in the evolution of the system state, plans of higher quality can be achieved by optimizing for the expected reward. In general, uncertainty in the model is represented by belief distributions over the system state and actions with probabilistic outcomes.

This section uses the notation of [Hem66]: Underlined symbols denote random variables whose realizations follow some probability distribution. The same symbol may occur underlined and non-underlined. Here, this can be the distinction between a random variable and a sampled realization of the random variable. Another case is the distinction between an action used as an identifier and the same action used to sample stochastic state transitions.

The system state following an initial state σ and an action with uncertainty *a* is in effect the realization of a random variable σ <sup>0</sup> ∼ σ <sup>0</sup> = *a*(σ). The probability (density) of the possible outcomes σ <sup>0</sup> ∈ Σ is P(σ <sup>0</sup> = σ 0 ). We allow the application of uncertain states with a belief distribution to actions The resulting state has a belief distribution P(*a*(σ) = σ 0 ) = R Σ [P(*a*(ν) = σ 0 ) P(σ = ν)] *d*ν. We do not consider the case where *a* is not applicable to the realization of the initial state σ.

By Equation 2.3, non-participating components *c /*∈ *C<sup>a</sup>* are not affected by the action execution. Actions with uncertainty are not exempt from this rule. All state transitions where a non-participating component changes its timed-indexed state via execution of an action *a* must have zero probability.

$$\forall \sigma \in \Sigma\_a, \forall \sigma' \in \Sigma, \forall c \notin C\_a: \sigma\_c \neq \sigma'\_c \Rightarrow \mathbb{P}(\underline{a}(\sigma) = \sigma') = 0 \tag{2.5}$$

If the state following the execution of an action with uncertainty is not immediately and fully knowable, we have to indirectly infer probabilities over the resulting state via incomplete or noisy observations. This is called the *partially observable* setting [KLC98]. With partial observability, the components *c* ∈ *C* each generate

²The most common reasons for delays and unforeseen events are, according to [VHL03], machine failure, urgent job arrival, job cancellation, due date change, shortage or delay in the arrival of material, change in job priority, rework or quality problems, over- or underestimation of process time, and operator absenteeism.

observations *o<sup>c</sup>* ∈ *Oc*. An action results in observations from the participating components *O<sup>a</sup>* = (×*c*∈*<sup>C</sup> Oc*). With a slight abuse of notation, the next state and the observations are both drawn by sampling the action (σ 0 *, oa*) ∼ *a*(σ). The distribution of observations conditionally depends on the action and the state transition P(*o<sup>a</sup>* |σ*, a,*σ 0 ).

Obviously, for every action *a*, the observations *o<sup>a</sup>* are conditionally independent of the state of components that do not participate in *a*. Otherwise, if an observation was conditionally dependent on the state of a non-participating component, then the observation could depend on the outcome of an action that occurs earlier in the action sequence but actually has a later starting time in simulation. This cannot be possible. Therefore, given the state of the participating components σ*a*, the resulting state of the participating components and the observations are conditionally independent of the non-participating components.

$$(\sigma'\_a, o\_a \mid \sigma\_a) \perp \sigma\_{C \backslash C\_a} \tag{2.6}$$

Now we extend the treatment of uncertainty to action sequences. We write Θ = ∪*a*∈*A*Θ*<sup>a</sup>* for the set of possible action parameters and *O* = ∪*a*∈*AO<sup>a</sup>* for the set of possible observations. The action-subscript for parameters and observations is dropped when the relation is clear from context.

**Definition 2.8.** *A history* h *is a sequence of episodes h k . Each episode consists of the selected action and action parameters, as well as the generated observation.*

$$h = \underbrace{a^1 \theta^1 o^1}\_{h^1} \underbrace{a^2 \theta^2 o^2}\_{h^2} \dots \underbrace{a^{|h|} \theta^{|h|} o^{|h|}}\_{h^{|h|}}$$

The set of possible histories *H* ⊂ (*A* × Θ × *O*) ∗ implies a tree-graph. The histories in *H* have finite length. Either because the scenario is *done* when a specific state is reached or by a cutoff at a maximum history depth. A history h can be appended with the next action *a*, parameters *θ* and observation o to form h <sup>0</sup> = h*aθ*o. In the case of a partially observable system, the current system state can only be inferred with uncertainty. Every history yields a probability distribution for the belief over the final system state that is conditioned on the

observations. Recall that *a*˜ denotes the already parameterized actions.

$$\sigma' \sim \underline{\hbar}(\sigma), \,\underline{\hbar}(\sigma) = \mathbb{P}\left( (a^{1\theta^1} \cdots a^{|\underline{h}|\cdot\theta^{|\underline{h}|}}) (\sigma) \, \middle| \, a^1, \dots, a^{|\underline{h}|} \right) \tag{2.7}$$

The computation of the belief distribution for the resulting state can be performed with iterative Bayes updates [Jay03] for the intermediary system state between actions. The state belief distribution following from an uncertain initial state is h(σ). In general, this computation cannot be performed with the available computational resources. The planning algorithms from the later chapters therefore rely on samples from forward simulations of the system only.

## **2.4 Reward and Policies**

So far, only the dynamics of actions and action sequences were discussed. Now we begin to express preferences between action sequences.

**Definition 2.9.** *Executing action a with parameters θ inducing a transition of the system state from* σ *to* σ <sup>0</sup> *generates a reward r*(σ*, a, θ,*σ 0 )*. The reward function is*

$$
\mathfrak{r}: \Sigma \times A \times \Theta \times \Sigma \to \mathbb{R} \,. \tag{2.8}
$$

With the reward function, each state defines a planning problem with the goal to maximize the future (expected) reward. This is known as the decision-theoretic planning problem [De 70; BDH99]. The following hierarchy of planning problems is distinguished in the literature [LaV06]:

The application of decision-theoretic planning for the (feedback) control of dynamical systems is discussed in [DW91; Ber+95; BG01]. See [LP12] for an application of MDP solvers to model manufacturing scenarios under uncertainty.

Algorithms for solving decision-theoretic planning problems belong to two distinct groups, online planning algorithms and policy constructing algorithms. Online planning algorithms are executed for the selection of the next action [Ros+08; SV10] at runtime. Policy constructing algorithms are executed ahead of time. They compute a fixed policy function that takes the current system state (or observed history) to the next action [KHL08].

**Definition 2.10.** *A policy is a mechanism to select the next action during runtime. In fully observable settings, policies are represented as functions π* : Σ → *A that map from the current system state to the next action. The policy function becomes π* : Σ → *A* × Θ *in settings with parameterized actions.*

*In POMDP, policies π* : *H* → *A or π* : *H* → *A* × Θ *map from the observed history to the next action. Internally, the POMDP policy may decide the next action and parameters based on a belief distribution over the current system state conditioned on the observed history* P(σ ∈ Σ | h ∈ *H*)*.*

For simplicity of the exposition, consider fully observable settings without parameterized actions. For a fixed policy *π* and a discount rate *γ* ∈ [0*,* 1), the value of an initial state σ 0 is the expected discounted reward.³

$$\mathfrak{v}^{\pi}(\sigma^{0}) = \mathbb{E}\left[\sum\_{k=0}^{\infty} \gamma^{k} \mathfrak{r}(\sigma^{k}, a^{k}, \sigma^{k+1}) \, \Big|\, a^{k} = \pi(\sigma^{k}), \, \sigma^{k+1} \sim a^{k}(\sigma^{k})\right] \tag{2.9}$$

As an alternative to the discount factor, we can limit the number of considered periods. The optimal value of a state v is the expected reward resulting from optimal action selection in every step. Based on Bellman's principle of optimality [Bel57], the V-value can be written as a recursive formula.

$$\mathfrak{v}(\sigma) = \max\_{a \in A} \underset{\sigma' \sim a(\sigma)}{\mathbb{E}} \left[ \mathfrak{r}(\sigma, a, \sigma') + \gamma \mathfrak{v}(\sigma') \right] \tag{2.10}$$

³The V-value and Q-value are longstanding terms in the literature. As both are scalar values we denote them as v and q in the mathematical notation.

The corresponding *Q*-value is the expected reward for selecting action *a* in state σ and optimal decision making thereafter.

$$\mathfrak{q}(\sigma, a) = \underset{\sigma' \sim a(\sigma)}{\mathbb{E}} \left[ \mathfrak{r}(\sigma, a, \sigma') + \gamma \mathfrak{v}(\sigma') \right] \tag{2.11}$$

Selecting a value according to max*a*∈*<sup>A</sup>* q(σ*, a*) is the optimal policy. But computing the optimal policy is generally computationally intractable. The remainder of this thesis is concerned with ways of rendering optimization and planning in the framework introduced in this chapter computationally feasible. This is achieved a) with tailored planning algorithms, b) so-called rollout policies that exploit known structure in the planning problem and c) the decomposition of the planning problem into a coupled set of smaller problems that are jointly optimized by cooperating agents.

# **3 Simulation-Based Planning for Concurrent Production Systems**

*Programming, or program planning, may be defined as the construction of a schedule of actions by means of which an economy, organisation, or other complex of activities may move from one defined state to another, or from a defined state toward some specifically defined objective.*

*Marshal K. Wood and George B. Dantzig [WD51]*

The model from Chapter 2 is generic and can be used to represent many types of systems. Very little constraints are imposed on the system dynamics that can be represented. This chapter develops an algorithm for sequential decision making that does not make additional limiting assumptions. But this richness with respect to possible system dynamics is a drawback when it comes to planning. Most planning algorithms impose a much more limited model structure they exploit to reduce the computational effort. The core postulate of this chapter is the following:

> *The same algorithm can be used for planning and runtime control on all levels of the control hierarchy – ranging from continuous dynamics of a physical system to global supplychain operations – and for both continuous and discrete production.*

This chapter develops a planning algorithm for the full model from Chapter 2 without additional assumptions. In the first two sections, two techniques are used to reduce the number of action sequences that are visited for planning with action

sequences in deterministic scenarios. In Section 3.1 the search tree is explicitly pruned by removing equivalent action sequences (this will have a precise definition). Section 3.2 further speeds up planning by implicitly pruning less promising parts of the search tree via Monte-Carlo Tree Search (MCTS). Section 3.3 extends planning to the full model with parametric actions and uncertainty with partial observability. In order to scale to larger scenarios, Section 3.4 develops a custom rollout policy that uses a relaxation of the planning problem to a Mixed-Integer Linear Program (MILP).

## **3.1 Tree Search with Backtracking**

The model from Chapter 2 can represent concurrency, i.e. parallelism and the synchronization of system components. A consequence of the model is that many action sequences are equivalent with respect to their overall preconditions and effects.

**Example 3.1.** Take the scenario from Example 2.2 with the modification that the robotic manipulator initially contains a product. The actions produce and package do not share a participating component as *C*produce = {mt} and *C*package = {manip}. It is easy to see that the sequence (produce package) is equivalent to the sequence (package produce) in terms of their preconditions and effects. In that sense, the two actions are independent from each other.

This notion of action sequence equivalence is made rigorous based on established results from Trace Theory (TT) [CF69; Maz77]. The theoretical prerequisites are summarized in the following. Full proofs can be found in [DM97a]. Some elementary definitions from Group Theory [DF04] are assumed as known.

### **Background: Trace Theory**

Previous nomenclature from Chapter 2 applies to TT by taking *A* as a set of letters of action-identifiers that are concatenated to words corresponding to action

sequences. Let *A* be a finite alphabet of letters and *A*<sup>∗</sup> the set of all finitelength words over *A*. The composition of letters forms a free monoid with generating set *A* and the empty word, denoted *ε*, as the unit element. Some letters commute and are said to be mutually independent. This is captured in the set of independence relations *Z* ⊆ *A* × *A*. Independence relations are symmetric (*a, b*) ∈ *Z* ⇒ (*b, a*) ∈ *Z* and irreflexive ∀*a* ∈ *A,*(*a, a*) ∈*/ Z*. We also write *a* ⊥*<sup>Z</sup> b* to denote independence between *a* and *b*. From *Z* follows an equivalence relation between words '*Z*: two words w and v from *A*<sup>∗</sup> are equivalent according to '*<sup>Z</sup>* if v is a permutation of w that can be reached by successive reordering of adjacent elements that commute according to *Z*. In algebraic terms, the trace monoid M(*A, Z*) is the quotient of *A*<sup>∗</sup> by the congruence '*Z*. Its elements, called traces, are pairwise disjoint subsets of *A*<sup>∗</sup> . Each trace contains the words that are mutually equivalent according to '*Z*. Let [w]*<sup>Z</sup>* denote the trace generated from the word w. We can now rephrase equivalency as w and v generating the same trace w '*<sup>Z</sup>* v ⇔ [w]*<sup>Z</sup>* = [v]*Z*. When it is clear from context which independence relations apply, we simply write [w] for the trace and w ' v to denote equivalence.

The independence relations define a partial order between some of the letters in a word. This partial order is a "must appear before" binary relation ≺.

$$w^k \prec w^l \Leftrightarrow k < l \land \neg (w^k \perp w^l) \tag{3.1}$$

The transformation between equivalent words w ' v is described by a permutation of the element indices *τ* . The partial order of the word elements is invariant to this permutation and therefore

$$\forall k, l \in \{1, \ldots, |w|\}, \; w^k \prec w^l \Rightarrow v^{\tau(k)} \prec v^{\tau(l)}.\tag{3.2}$$

The partial order yields the same *dependence graph* and *Hasse diagram* for every word from the same trace.¹ The dependence graph G = (*V, E*)is a directed acyclic graph with nodes *V* = {*w* 1 *, . . . , w*<sup>|</sup>w<sup>|</sup>} and edges *E* = {(*w, v*) ∈ *V* 2 : *w* ≺ *v*}.

¹Graphs being equal refers to the existence of an isomorphism between the graphs that takes node labels into account.

**Figure 3.1:** The Dependence Graph (left) and Hasse Diagram (right) for the trace [*abcaed*] according to the trace monoid M(*A, Z*) over the alphabet *A* = {*a, b, c, d, e*} and independencies *Z* = {(*a, c*)*,*(*c, a*)*,*(*b, d*)*,*(*d, b*)*,*(*d, e*)*,*(*e, d*)}.

The Hasse diagram G <sup>0</sup> = (*V, E*<sup>0</sup> ) is the dependence graph with redundancies removed *E*<sup>0</sup> = {(*w, v*) ∈ *E* : @*u, w* ≺ *u* ≺ *v*}. See Figure 3.1 for an example.

Besides their dependence graph and Hasse diagram, traces are uniquely represented by a word in a normal form. That is, exactly one word from every trace is in the normal form. There exist several normal forms with this property. We consider the Lexicographical Normal Form (LNF): Assume a total ordering of the letters *A* and a resulting lexicographical ordering over words. Longer words have a larger order than smaller words regardless of their elements. For two words w*,* v of the same length, w *<* v if *w <sup>k</sup> < v<sup>k</sup>* for the smallest index *k* where *w <sup>k</sup>* 6= *v k* . A word is in LNF if it has the smallest lexicographical order among all words of the trace.

#### **Tree-Search with Trace-Based Pruning**

With the prerequisites in place, we now take a look back at Chapter 2. Remember that action sequences are operators on the system state with defined preconditions (the operator domain) and effects.

**Definition 3.1.** *Two action sequences* w *and* v *are equivalent* w ' v *if they have the same domain of initial states* Σ<sup>w</sup> = Σ<sup>v</sup> *and yield identical results* ∀σ ∈ Σw*,* w(σ) = v(σ)*.*

The independence of actions is defined based on the possibility to commute them if they occur in adjacent positions in any action sequence.

**Definition 3.2.** *Two actions a and b are independent a* ⊥ *b if commuting them yields an equivalent sequence* ∀w*,* v ∈ *A*<sup>∗</sup> *,* w*ab*v ' w*ba*v ⇔ *a* ⊥ *b.*

It is generally intractable to show the equivalence of two action sequences by enumerating the set of valid initial system states and the effect of the action sequence on those states. For some action pairs however, independence can be shown without evaluating the preconditions and effects:

**Proposition 3.3.** *Any two actions a, b that do not share a participating component are independent.*

$$C\_a \cap C\_b = \mathcal{Q} \Rightarrow a \perp b \tag{3.3}$$

*Proof.* For all action sequences w*,* v there must be w*ab*v equivalent to w*ba*v. This follows directly from the operational semantics of actions defined in Equation 2.3.

**Example 3.2.** Figure 3.2 shows which sequences from Example 2.3 are equivalent and can be pruned. Only the sequences up to a length of five are shown. The initial lexicographic ordering of the actions for the LNF is produce *<* put *<* take *<* package.

**Figure 3.2:** Sequence tree from Example 2.3 before (left) and after the pruning of equivalent sequences (right).

**Algorithm 1** Is the action sequence w*a* in LNF given that w is in LNF?

Some sequences may additionally be equivalent due to some lucky alignment of the action effects. Such additional equivalencies are not considered in this text. Now that the independence of any two actions is easily computed, we use this information to speed up the search for the best action sequence. Only one sequence from each trace is considered. The sequences from the same trace are equivalent also with respect to the reward that they generate. To consider only one sequence per trace, we evaluate only the action sequences in LNF.

**Proposition 3.4.** *If a word* w *is not in LNF, no word* wu *starting with the prefix* w *is in LNF.*

*Proof.* If w is not in LNF, then a word v ' w must exist so that v *<* w. For all u ∈ *A*<sup>∗</sup> there is vu ' wu and vu *<* wu.

From Proposition 3.4, we know that a sequence w*a* can only be in LNF if w is in LNF. During a depth-first traversal of the search tree, entire subtrees can thus be pruned away. If the action sequence w leading to the current position in the search tree is known to be in LNF, the candidate sequence w*a* can be tested for LNF very fast. The test is performed by Algorithm 1. It has a worst-case runtime that is linear in the sequence length. In practice the algorithm performs much faster as a breaking condition is usually found within the first few elements of the sequence.

**Proposition 3.5.** *If the action sequence* w *is in LNF, then for all a* ∈ *A Algorithm 1 returns true if and only if the sequence* w*a is also in LNF.*

*Proof.* For a sequence v, denote with *Z*(v*, k*) = {*l* : *l < k* ∧ ∀*m* ∈ {*l, . . . , k* − 1}*, v<sup>k</sup>* ⊥ *v <sup>m</sup>*} the indices of the contiguous elements before *v k* that all commute with *v k* . In the following we show that an action-sequence v is in Lexicographical Normal Form (LNF) if and only if

$$\forall k \in \{1, \ldots, |\mathbf{v}|\}, \forall l \in Z(\mathbf{v}, k), \ v^l < v^k. \tag{3.4}$$

First, we show that (3.4) must be satisfied for every sequence v in LNF. Assume v is in LNF and does not hold condition (3.4). Then there exists a combination of indices *l < k* where *v <sup>l</sup> > v<sup>k</sup>* and ∀*m* ∈ {*l, k* − 1}*, v<sup>m</sup>* ' *v k* . The decomposition v = u*v <sup>l</sup>*q*v <sup>k</sup>*z can be rearranged as y = u*v kv <sup>l</sup>*qz with possibly empty subsequence u*,* q and z. This contradicts v being in LNF since v ' y and y *<* v.

Second, we show that condition (3.4) is sufficient for v to be in LNF. Assume v satisfies (3.4) and is not in LNF. Then there must exist an equivalent sequence y ' v in LNF. Equivalent sequences are permutations and |y| = |v|. So for y to have a smaller order, there must be an index *l* where y and v first differ y :*l*−<sup>1</sup> = v :*l*−1 and *v <sup>l</sup> < y<sup>l</sup>* . Due to y ' v, action *v <sup>l</sup>* must appear in v at an index greater than *l* and v can be decomposed as v = u*v <sup>l</sup>*q*y <sup>l</sup>*z, with possibly empty subsequences u*,* q and z, where *y l* commutes with *v l* and all elements of q. This contradicts (3.4).

Since Equation 3.4 only refers to the preceding elements in the action sequence and w is known to be in LNF, it suffices to test (3.4) for the new element.

**Proposition 3.6.** *For any action a and LNF action sequence* w*, Algorithm 2 returns an action sequence that is in LNF and equivalent to* w*a.*

*Proof.* When the condition in line 3 of Algorithm 2 is not true, then *a* commutes with *w k* and *a < w<sup>k</sup>* . Therefore, once the condition in line 3 is true, w:*<sup>k</sup>a* is in LNF and the LNF-order of the sequence cannot be improved by moving an element from w*<sup>k</sup>*+1: before *a*. Since w is in LNF, no permutation of the sub-sequence w*<sup>k</sup>*+1: has reduced LNF-order. Therefore w:*<sup>k</sup>a*w*<sup>k</sup>*+1: is in LNF.

A trivial approach to search for good sequences is to cast the search space as a tree structure and to enumerate the solutions at the leaf nodes of the tree. Every node in the tree represents a (partial) action sequence. Depth-first search completely enumerates the solutions at the leaf nodes without having to hold the entire tree in memory. Algorithm 3, called DFS for depth-first search, walks over the tree in recursive fashion. The algorithm backtracks to a prior partial action )

**Algorithm 3** Depth-First Search

```
1: v
   max ← −∞, wmax ←
2: procedure DFS(ω, v, w)
3: if |w| = k
             max then
4: if v > vmax then
5: v
            max ← v
6: wmax ← w
7: return
8: for a ∈ A : ω ∈ Σa do
9: ω
          0 ← a(ω)
10: v
         0 ← v + r(ω, a, ω
                       0
11: DFS(ω
              0
              , v0
                 , wa)
12: DFS(σ, 0, )
13: return wmax
```
#### **Algorithm 4** Pruned Depth-First Search

```
1: v
   max ← −∞, wmax ←
2: procedure DFSᴸᴺF(ω, v, w)
3: if |w| = k
             max then
4: if v > vmax then
5: v
            max ← v
6: wmax ← w
7: return
8: for a ∈ A : ω ∈ Σa ∧ TestLNF(w, a) do
9: ω
          0 ← a(ω)
10: v
         0 ← v + r(ω, a, ω
                        0
                        )
11: DFSᴸᴺF(ω
                 0
                 , v0
                    , wa)
12: DFSᴸᴺF(σ, 0, )
13: return wmax
```
**Algorithm 5** Branch & Bound

```
1: v
   max ← −∞, wmax ←
2: procedure B&B(ω, v, w)
3: if |w| = k
             max then
4: if v > vmax then
5: v
            max ← v
6: wmax ← w
7: return
8: for a ∈ A : ω ∈ Σa do
9: ω
          0 ← a(ω)
10: v
         0 ← v + r(ω, a, ω
                       0
                        )
11: if v
           0 + β(ω
                 0
                  , wa) > vmax then
12: B&B(ω
                 0
                  , v0
                    , wa)
13: B&B(σ, 0, )
14: return wmax
```
#### **Algorithm 6** Pruned Branch & Bound

```
1: v
   max ← −∞, wmax ←
2: procedure B&BᴸᴺF(ω, v, w)
3: if |w| = k
              max then
4: if v > vmax then
5: v
            max ← v
6: wmax ← w
7: return
8: for a ∈ A : ω ∈ Σa ∧ TestLNF(w, a) do
9: ω
          0 ← a(ω)
10: v
         0 ← v + r(ω, a, ω
                        0
                        )
11: if v
           0 + β(ω
                 0
                  , wa) > vmax then
12: B&BᴸᴺF(ω
                    0
                    , v0
                       , wa)
13: B&BᴸᴺF(σ, 0, )
14: return wmax
```

sequence once all solutions starting at the current position in the tree have been enumerated. The tree of depth *k* max contains up to P*<sup>k</sup>* max *<sup>k</sup>*=0 |*A*| *<sup>k</sup>* nodes and |*A*| *k* max leaf nodes. Enumerating all possible combinations is intractable for all but the most trivial scenarios. Usually the branching factor of the tree is reduced because not all actions from *A* are applicable in the intermediary system states. But the problem of combinatorial explosion of the search space remains. Proposition 3.4 can be used to test at every node whether the entire subtree starting at the node is not in LNF and can be skipped (pruned). Since every trace contains exactly one sequence in LNF, it is ensured that traversing the tree without the pruned branches still visits every trace once. Algorithm 4, called DFSᴸᴺF, extends DFS with the pruning of sequences that are not in LNF. Note that only line 10 of Algorithm 4 has changed compared to DFS.

Another technique to speed up tree search is Branch & Bound [Lit+63]. Branch & Bound is a well-known technique for combinatorial optimization in the presence of so-called admissible heuristics [Pea84]. A heuristic is called admissible if it overestimates the performance that can still result from a partial solution. Here, we express the admissible heuristics as a function *β*(σ*,* w) with the current system state σ and the current partial action sequence w as arguments. If the admissible heuristics show that the best performance starting from the current partial solution is worse than the best solution that was already encountered, then the entire subtree behind the current partial solution can be pruned. Branch & Bound is implemented in Algorithm 5. In Algorithm 6, Branch & Bound is combined with trace-based pruning. Again, only line 10 is changed to prune out sequences that are not in LNF. The commonality between DFS and Branch & Bound is the use of backtracking to return to a previous partial solution either when all nodes in a subtree have been visited or when the solutions that remain in the subtree are known to be suboptimal.

#### **Evaluation**

The effect of the pruning techniques for tree search methods is evaluated based on the Jobshop Scheduling Problem (JSP) [Pin08]. Scheduling is one of the most important planning problems on the shopfloor. It is also computationally challenging. The jobs *j* ∈ *J* have to be processed on the machines *M*. Every job consists of operations *o* ∈ *O<sup>j</sup>* = {1*,* 2*, . . . ,* |*M*|}. The operations need to be processed in-order for every job and are each assigned to a specific machine *mj,o* ∈ *M*. The operation duration is *dj,o*. The goal is to find a schedule that assigns operations to machines in order to minimize the finishing time of the last job. See Table 3.1c for a benchmark JSP from the literature. There are (|*J*|!)|*M*<sup>|</sup> possible schedules for a JSP with |*J*| jobs and |*M*| machines where every jobs needs to visit every machine once [JM98]. So a JSP with |*J*| = 20 and |*M*| = 10 has 7*.*2651 × 10<sup>183</sup> possible solutions. Compare this to the mass of the observable universe currently estimated to the equivalent of 10<sup>80</sup> hydrogen atoms. Even though the number of possible solutions is huge, current solution techniques can compute near-optimal solutions for JSP with hundreds of machines and jobs. Many benchmark scheduling problems have even been solved with certified optimaltiy. However, the worst-case analysis of the JSP shows that it is NP-hard for instances where |*J*| ≥ 3 and |*M*| ≥ 3 [GJS76]. Therefore, unless P = NP, all algorithms for solving the JSP are either heuristic or require a long running time for some (synthetic) JSP problems.

We cast the JSP in the model from Chapter 2 in terms of actions with preconditions and effects. Every machine and every job in a JSP is represented as a component in the system model *C* = *C*machines ∪ *C*jobs. The operations for each job are represented as actions where both one machine and one job participate. The time-indexed state of the machine components is trivial. It only consists of the time when the machine is available next. The state of the job components is the number of operations that were already performed for the job. The precondition of every action is that the associated operation is next in line for the job. The action effect simply increases the number of finished operations for the job by one. The duration of each action (operation in the JSP) is deterministic. The cost generated by an action is the increase in the maximum timestamp of the components. So

**(c)** The 6 × 6 benchmark JSP ft06 from [FT63].

**Table 3.2:** Example JSP problems.

the cost of the entire sequence is the makespan of the solution.

$$\mathfrak{x}(\sigma, a, \sigma') = \left(\max\_{c \in C} t\_c\right) - \left(\max\_{c \in C} t'\_c\right) \tag{3.5}$$

Here, *t<sup>c</sup>* refers to the simulation time of component *c* in the system state σ and *t* 0 *c* refers to the simulation time of the component in the system state σ 0 .

Given a partial schedule as an action sequence w, the following trivial heuristic computes an upper bound for the reward that can be generated following the execution of w. For every component (representing a job or a machine), the remaining operations are summed up and added to its current time.

$$\beta^{JSP}(\sigma, w) = -\max\_{c \in C} \left( t\_c + \sum\_{\substack{a \in A : a \notin w, \\ c \in C\_a}} d\_a \right) \tag{3.6}$$

The ft06 benchmark problem from Table 3.1c illustrates the importance of optimisation. An unoptimized solution can require more than 150 seconds to complete all jobs. The optimal solution requires just 55s. Finding the optimal solution for larger problems is not possible without the help of computers. In addition to the ft06 benchmark, two additional minimal problems are considered. See Table 3.2 for their exact definition. The minimal examples are too small to be of any practical value. But they give an indication of how fast the problem complexity increases. A sequence of 9 unique elements has 9! = 362*,* 880 permutations. The minimal 3 × 3 JSP depicted in Table 3.2 also defines 9 actions. But the full scenario tree for the JSP has only 1*,*680 leaf nodes corresponding to valid schedules. By additionally pruning branches that are not in the LNF only 63 leaf nodes remain. Using the makespan lower bound from Equation 3.6, nodes can be pruned where the lower bound is equal or worse to the best solution encountered so far. With the Branch & Bound pruning, only 3 leaf nodes are actually visited. However, more nodes were visited overall compared to trace-based pruning for full depth-first search. This indicates that mostly branches close to the leaves were pruned with Branch & Bound only. By combining trace pruning with Branch & Bound, only 65 nodes are visited overall. This is a decrease in the search running time by a factor of more than 50 compared to Branch & Bound search without pruning.

On the 4 × 4 JSP depicted in Table 3.2 an even smaller fraction of the scenariotree is visited by combining trace-based and makespan-based pruning. However, the number of visited nodes still grows fast with the size of the JSP. On the ft06 6×6 JSP, 682,508 nodes are expanded with both bound and trace pruning enabled. Our implementation is capable of visiting over 1,000,000 nodes per second during Branch & Bound search (on a Lenovo T480 laptop computer). However, with only either trace-based pruning or Branch & Bound enabled, we could not solve ft06 in over a day of computation. A 10 × 10 JSP could not be solved within several days of compute time even with both trace-based and bound-based are activated.

#### 3 Simulation-Based Planning for Concurrent Production Systems

**(a)** Visited nodes for the minimal 3 × 3 JSP

**(b)** Visited nodes for the minimal 4 × 4 JSP

**(c)** Visited nodes for the ft06 6 × 6 JSP

**Table 3.3:** Benchmarking of pruning techniques for DFS and Branch & Bound. Questionmarks indicate that an optimization did not terminate after two weeks of computation. So the numbers are outstanding.

In summary, trace-based pruning reduces the number of visited action sequences to a small fraction. For the 4 × 4 JSP, the number of action sequences (visited leaf nodes) was reduced by a factor of more than 5,600. The improvements are comparatively increasing with the length of the action sequences.

We have shown that LNF pruning leads to a dramatic reduction of the search space. The reduction is already on many orders of magnitude for the comparably small benchmark problems that were considered. However, even with tracebased pruning, naive tree search does not scale up to scenarios of relevant size. In practice, genetic algorithms are often used to solve combinatorial problems. Modern heuristic solvers find good solutions for JSP with thousands of jobs [Dim15]. But these solvers exploit the specific structure of the JSP. Handling of complex action preconditions is near-impossible for genetic algorithms. The genetic crossing and mutation of two partial solutions will almost always lead to an infeasible action sequence where the preconditions of an action are not met. Repairing an infeasible plan is usually quite difficult. For example if a single product is missing to finish an order of many hundred products.

## **3.2 Planning for Discrete Action Sequences**

The algorithms from Section 3.1 use backtracking to return to previous partial solutions. This requires either that the system state is fully known so that it can be stored in a computer (to return to previous states) or that action sequences can be deterministically repeated. For experiments carried out in the physical world, neither of these assumptions is true. In this section, we further speed up planning by "implicitly pruning" less promising branches with Monte-Carlo Tree Search (MCTS). Many recent breakthroughs in Artificial Intelligence were made possible by MCTS. This includes the AlphaGo system [Sil+17] which is able to play the game of Go with superhuman performance. In short, MCTS enables planning in scenarios with many combinatorial variations where backtracking tree search is not able to cover any significant portion of the search space. In addition, MCTS uses only forward-simulation without backtracking. This reduces the coupling between the planning algorithm and the simulator to the point where real-world simulations could be used to generate samples for the planning algorithm.

#### **Background: Monte-Carlo Tree Search**

In MCTS [Bro+12; Mun+14], a scenario tree is explored by iterative playouts. A playout is essentially one run of the scenario with sequential decision making in a series of steps. Every playout starts at the root of the sequence tree and evolves by "forward simulation". In contrast to Branch & Bound, there is no backtracking to previous states within a playout. Historically, MCTS evolved from research on the multi-armed bandit problem: The multi-armed bandit problem is an idealized version of a slot machine. The following short exposition follows [Mun+14] and adapts the notation to the conventions of this text. Consider a multi-armed bandit in a casino where they player can choose a different arm during *k* max rounds (this will be extended to sequential decision making later on). The reward of the different arms *a* ∈ *A* is random according to the distribution P*<sup>a</sup>* with a support on [0*,* 1]. But the player is initially unaware of the reward distribution for every arm. In each round *k*, the player chooses a bandit *a <sup>k</sup>* ∈ *A* and collects a reward *r <sup>k</sup>* ∼ P*a<sup>k</sup>* . The rewards generated at the other arm are not observed. The goal is to find a strategy that maximizes the expected sum of payouts during the *k* rounds. This leads to the so-called *Exploration/Exploitation Tradeoff* : The player could select a bandit with high expected reward. But this might get him stuck at a bandit with suboptimal expected reward. So the player wants to explore other options without loosing too much in the process. The expected reward of the different arms *a* ∈ *A* is *µ<sup>a</sup>* = E[P*a*]. The expected reward of the best arm (there can be several best arms) is *µ* <sup>∗</sup> = max*a*∈*<sup>A</sup> µa*. Possible strategies for repeated play are analyzed with respect to their *expected regret*. The cumulative regret after *k* rounds compares the actual reward *r* with the expected reward of choosing the arm with highest expected reward every time.

$$b\_k = k\mu^\* - \sum\_{l=1}^k r^l$$

The expected cumulative regret is therefore

$$\mathbb{E}[b\_k] = k\mu^\* - \sum\_{l=1}^k \mathbb{E}[r^l] = \sum\_{a \in A} \mathbb{E}[n\_a(k)](\mu^\* - \mu\_a)\,,$$

where *na*(*k*) denotes the number of draws from *a* after the first *k* rounds.

An important tool for bounding the expected cumulative regret is the Chernoff-Hoeffding inequality. Let (*y i* )*<sup>i</sup>*=1*,...,n* i.i.d. samples of a probability distribution with support [0*,* 1] and mean *µ*. The empirical mean estimator *µ*ˆ = 1 *n* P*<sup>n</sup> <sup>i</sup>*=1 *y i* is a function of the *y i* and therefore a random variable as well. The Chernoff-Hoeffding inequality gives probability bounds for the distance between the true and the estimated mean.

$$\mathbb{P}\left(\underline{\hat{\mu}}-\mu \geq \epsilon\right) \leq e^{-2n\epsilon^{2}} \quad \text{and} \quad \mathbb{P}\left(\underline{\hat{\mu}}-\mu \leq -\epsilon\right) \leq e^{-2n\epsilon^{2}}$$

The bound is independent of the underlying distribution of the *y<sup>i</sup>* . So it can be applied even when the distribution of the *y i* is not known!

There exists a variety of approaches for selecting the next draw in the multiarmed bandit setting. Auer et al. [ACF02] proposed the so-called Upper Confidence Bound (UCB) algorithm. Let *µ*ˆ*a,k* denote the mean return of arm *a* sampled during the first *k* rounds. UCB always selects the next arm as follows:

$$a^k = \operatorname\*{arg\,max}\_{b \in A} \hat{\mu}\_{b, k-1} + \sqrt{\frac{3 \log k}{2n\_b(k-1)}}$$

Selecting the best arm according to the UCB gives an upper bound for the expected cumulative regret by application of the Chernoff-Hoeffding inequality.

$$\mathbb{E}[b^k] \le 6 \sum\_{\substack{a \in A, \\ \mu^\* > \mu\_a}} \frac{\log k}{\mu^\* - \mu\_a} + |A|(\frac{\pi^2}{3} + 1)\,,$$

The expected cumulative regret grows at most logarithmically in *k*.

The algorithm of Auer et al. [ACF02] selects actions from a flat list of options. Monte-Carlo Tree Search (MCTS) is a family of algorithms that uses the same principle of iterative exploration for planning in sequential decision settings. A popular variation of MCTS, called Upper Confidence on Trees (UCT) [KS06], uses the UCB decision rule for every action choice during a playout. See Figure 3.3 for an overview. Algorithm 7 shows UCT (but with some modifications compared to [KS06] that will be explained in the following). MCTS can be seen

**Figure 3.3:** Outline of a Monte-Carlo Tree Search. Reproduced from [Cha+08].

as "simulation-based" planning, as no backtracking is used. Instead, the results from the last simulation are incorporated into statistics about the expected reward for the possible next actions at the current position in the scenario tree. The updated empirical reward statistics is used for decision-making in the following playouts.

The assumption that the reward in every step is in the range [0*,* 1] is not required for UCT to converge to optimal decisions in the limit. On the downside, there is currently no good characterization of the convergence speed achieved by UCT that doesn't make strong assumptions on the underlying scenario. It is only known that for a large enough number of plays, every branch of the scenario tree is visited. Since only a small subset of the tree is explored in practice before the algorithm is terminated, MCTS struggles with scenarios where the reward is "unevenly distributed". That is, if a big reward can be found after a long sequence of actions where small changes to the action sequence lead to much worse results, it is then not very likely that the reward is encountered at all and MCTS will not assign a correct value estimation for the choices. A workaround for this is Reward Shaping [NHR99] where the reward is distributed such that it may be encountered early on in the action sequences.

UCT may require tuning of the parameter *α* which regulates the importance of the upper confidence bias. Setting *α* corresponds to making a choice for the Exploration/Exploitation tradeoff between the two extremes "always explore" and

"always exploit". Note that UCB was developed for action selection in stochastic scenarios with unknown reward distributions for the different actions. Here we apply the same principle to tree-search for deterministic scenarios. There is an informal argument for this. Many rollout policies are stochastic, similar to the uniform sampling policy described above. So the rollout takes samples from the reward distribution that is implied by applying the current rollout policy to the subsequent steps. When nodes are visited several times, the underlying distribution for the policy changes. So past experience does not necessarily match the samples taken later on. Nevertheless, the UCB-based selection rule shows good performance in practice.

### **Monte-Carlo Tree Search for Discrete Action Sequences**

We introduce a modification of UCT called UCTᴸᴺF. See Algorithm 8 for details. UCTᴸᴺF is inspired from UCT but deviates in a few key aspects. First, the original UCT algorithm gradually builds up a tree structure where the nodes correspond to partial action sequences. But only one new node is added with every play. Once the edge of the current tree is reached, the remaining steps are "rolled out" and only the sum of rewards for the rollout is considered. UCTᴸᴺF does not use rollouts and adds nodes for all selected actions in the sequence. Second, the original UCT stores a statistics about the average reward that was achieved after selecting a node (action). Instead, we perform a maximization in every step of the Update procedure. Therefore, the value estimation of every node is the maximum reward that can be achieved by following the best known sequence of actions after selecting the node (action). Third, we do not only update a single sequence for every play. Instead, when we arrive at a sequence that is not in LNF, then we permute the sequence to the equivalent LNF sequence. But we store the original sequence and the LNF sequence and run the Update procedure on both of them. All three modifications to the original UCT algorithm are interrelated. This is explained next.

The easiest way to restrict the search to sequences in LNF only is to simply remove actions leading to non-LNF sequences at the current node in the search tree. But when choosing one action at a time, this often leads to dead-ends where no action can be chosen without breaking the LNF constraint. Since we do not want to backtrack to previous system states during a playout, virtually no playout would complete. Therefore we use the LNF algorithm to repair the action sequence as we go along. However, simply permuting the sequence to LNF after every action selection is not compatible with the UCT approach either: Soon, some nodes are always selected because they have a very low visit count. But as the sequence is permuted to LNF, this nodes continue to not have their visit count *n* increased. Suppose an example with three actions *a < b < c* that all commute and where every action can be chosen only once in every sequence. The sequence starting with the action *c* would always become the LNF sequences *ac* or *bc* after the second action choice. So the node for the first action *c* is nevery updated. But the UCT algorithm will always choose *c* first as long as the counter *n*[*c*] is not increased. This is solved by updating all sequences that were encountered in the current playout, regardless of whether they are in LNF or not. The counter *z* for the node updates within the current playout ensures that no node is updated twice for one playout. Note that we use the algorithm LNF in a slightly different fashion. Giving the action sequence w and the reward history r as input, both

are permuted to yield an LNF action sequence where the index of actions still matches the index of the corresponding reward.

#### **Evaluation**

We evaluate the algorithms UCT and UCTᴸᴺF on the benchmark Jobshop Scheduling Problems (JSP): The ft06 benchmark with 6 jobs on 6 machines from [FT63] and the abz5 benchmark with 10 jobs on 10 machines from [ABZ88]. However, we change the classical JSP in one important aspect: Instead of optimizing the makespan, the time when the last job finishes, we target the sum of the finishing times (the tardiness if the job is immediately due) for all jobs. The reason for the change is to ensure that the immediate reward of actions depends on a small number of participating components *C<sup>a</sup>* only. When optimizing the makespan, it is of course possible to return the negative increase of the maximum simulation time for all components. But then the actions would have all components *C* as participants. We prefer to have only a small subset of the components participating in each action, so that the trace-based equivalence of action sequences as defined in Section 3.1 can be used to prune the sequence tree. The rollout policy *π* is set to select randomly among the possible choices with a uniform distribution.

As can be seen in Figure 3.4, MCTS makes big improvements to the best-known solution early on. Depending on the exploration parameter *α*, the algorithm then "converges" towards a reward value that is no longer improved upon even with very long running times. This is not a true convergence however, since we know that all branches of the search tree are explored eventually with the UCB rule. The more complex abz5 example shows an interesting phenomenon where higher *α* lead to better results. But it takes longer until "convergence" is reached. This insight makes sense with regards to the UCB action selection. Once all actions at a given node have been explored several times, the UCB rule will predominantly choose actions with a known high reward. So the observed convergence is explained by a shift from exploration to exploitation in the action selections. Pruning sequences that are not in LNF speeds up the convergence to some final best reward. This was expected as the pruning effectively reduces the size of the search tree.

#### 3 Simulation-Based Planning for Concurrent Production Systems

**(a)** Benchmark of UCT with the JSP example ft06.

**(c)** Benchmark of UCT with the JSP example abz5.

**(b)** Benchmark of UCTᴸᴺF with the JSP example ft06.

**(d)** Benchmark of UCTᴸᴺF with the JSP example abz5.

**Figure 3.4:** JSP benchmarks for Monte-Carlo Tree Search. Every benchmark was run 10 times. The lines give the average empirical reward. The standard error is indicated by tick marks and the variance is shown in a lighter shade.

#### **Algorithm 8** UCT for Deterministic Actions with Trace-Based Pruning

```
1: procedure UCTᴸᴺF(σ
                   0
                    )
2: n[ · ] ← 0, q[ · ] ← 0
3: while enough time do
4: Y ← PlayᴸᴺF(σ
                     0
                      )
5: UpdateᴸᴺF(Y )
6: return arg max
             a∈A
                  q[a]
1: procedure UpdateᴸᴺF(Y )
2: z[ · ] ← 0
3: for (w, r) ∈ Y do
4: for k = |w| . . . , 1 do
5: g ← w:k
6: if z[g] > 0 then
7: break
8: n[g] ← n[g] + 1, z[g] ← 1
9: q[g] ← r
                  k + max
                      a∈A,
                     n[ga]>0
                           q[LNF(ga)]
1: procedure PlayᴸᴺF(σ)
2: w ← ε, r ← ε, Y ← ∅
3: while ¬done(σ) do
4: B ← {b ∈ A : n[wb] = 0}
5: if B 6= ∅ then
6: a ← πB(σ)
7: else
8: a ← arg max
               b∈A, σ∈Σb
                       h
                        q[LNF(wb)] + α
                                     q log n[w]+1
                                       n[LNF(wb)] i
9: (σ, v) ← a(σ)
10: Y ← Y ∪ {(wa, rv), LNF(wa, rv)}
11: (w, r) ← LNF(wa, rv)
12: return Y
```
# **3.3 Planning with Uncertainty and Continuous Action Parameters**

As defined in Section 2.2, parameterized actions a take parameters from the set Θ*a*. But the algorithms from the previous sections 3.1 and 3.2 can only handle discrete choices of deterministic actions. If the set of parameters is finite, we can simply extend the tree-search techniques to search over the joint space of actions and their parameters. This is however not possible if the action parameters are uncountable. For example if a parameter is defined on a continuous domain. Furthermore, while the UCB rule for bandit problems is defined for decision making in stochastic unknown environments, UCT assumes deterministic actions. This section deals with the extension of the simulation-based planning approach to actions *a* that can be both stochastic and have parameters on a continuous parameter space.

## **Background: MCTS under Uncertainty**

For MDP, where the outcome depends not only on decision-making, but also on the response from the stochastic system model, it has been known for some time that a randomly sampled subset of the scenario tree that covers only a vanishing fraction of the full scenario is enough to compute near-optimal actions from any state [KMN02]. If is the admissible error for the estimation of the V-value of the current state (see Section 2.4), then the number of sample playouts grows exponentially in . It does however not depend on the size and complexity of the state representation. With the advent of MCTS, sampling based methods have also been used for sequential decision-making under uncertainty [KS06].

The POMCP algorithm (Partially Observable Monte-Carlo Planning) extends MCTS to scenarios with partial observability (POMDP) [SV10]. The authors [SV10] write with regards to the performance of their invention POMCP:

[On a benchmark problem], POMCP achieved the same performance with 4 seconds of online computation to the state-of-the-art solver SARSOP with 1000 seconds of offline computation.

While the current system state σ is known to the simulator used to sample state transitions, the decision-making only relies on the observations resulting from the

actions and the reward statistics that were collected prior. See Algorithm 9 for details. The search tree now not only consists of nodes representing actions, but is a bipartite graph of actions and the resulting observations. The rollout policy *π* takes the full history of actions and observations as input. A trivial rollout policy is to sample uniformly from the previously unexplored actions. This policy is a popular choice as it does not depend on prior knowledge and assumptions about the scenario and the state representation.

## **Background: Optimistic Optimization**

Optimization of functions on a continuous domain has a long history. Very efficient solvers exist today for the optimization of convex functions where gradient information is available [BV04]. Gradient-free global optimization of non-convex functions remains challenging. Classical solution techniques are the Dividing Rectangles (DIRECT) algorithm [JPS93] as well as heuristic genetic algorithms

**Figure 3.5:** Example for Optimistic Optimization on the unit cube. The unit cube is split into cells based on a hierarchical partitioning of the domain. Every cell has an index (~*, i*) that also represents a node in the partitioning tree: ~ is the height of the tree at that node (the number of splits) and *i* is the index within the set of cells of height ~. The initial node (0*,* 1) represents the entire domain, here depicted as the unit cube. Leaf nodes are marked gray in the partitioning tree.

[SP97]. While the latter provides no lower bounds for performance, the convergence speed of DIRECT depends on an upper bound of the function's gradient known as the Lipschitz bound. In recent years, the ideas for optimization in discrete bandit settings and MCTS have been translated to optimization of continuous functions. The resulting techniques are known as Optimistic Optimization (OO). See [Mun+14] for a comprehensive account.

Suppose we want to maximize some function *f* : *X* → R over the domain *X*.

$$x^\* = \underset{x \in X}{\text{arg}\max} \, f(x) \tag{3.7}$$

If we can make assumptions of convexity, then a range of established methods and commercial tools can be used to solve the optimization problem [BV04]. OO does not require an assumption of convexity. The basic idea is to iteratively dissect *X* into disjoint cells of decreasing size that cover the entire function domain (Figure 3.5a). The hierarchy of cells forms a tree structure (Figure 3.5b). An "optimistic" upper bound is computed for every cell. This upper bound guides the continuing selection and splitting of the cells. The breakthrough of recent work is to not require a Lipschitz bound for *f* to compute the upper bound. Instead, the function is assumed to be smooth around the optimizer with respect to a semi-metric.² This is a much weaker assumption.

**Algorithm 10** Deterministic Optimistic Optimization (DOO) for the optimization of an unknown function *f* : *X* → R. This formulation is for the special case where the search domain is the *n*-dimensional unit cube *X* = [0*,* 1]*<sup>n</sup>* and cells are split into three children.

**Deterministic Optimistic Optimization** The first considered OO algorithm is Deterministic Optimistic Optimization (DOO). Algorithm 10 shows a simplified version of DOO. It assumes the domain of *f* is the *n*-dimensional unit cube and cells are split into three children. In every iteration, the cell with the highest upper bound is selected and split into 3 children. The cells are denoted as (~*, i*) where

²A semi-metric has *`*(x*,* y) = *`*(y*,* x) and *`*(x*,* y) = 0 ⇒ x = y. Different from a regular norm, the triangle inequality (a consequence of the Cauchy-Schwarz inequality) is not required to hold.

~ is the depth of the tree and *i* is an index for the cells at the same depth.³ The set *L* contains the leaf nodes of the current search tree (cf. Figure 3.5b). The midpoint of the visited cells (~*, i*) is stored as x[~*, i*]. The upper bound of each cell is computed from an evaluation at the midpoint and the bias *δ*(~) that depends on the depth-position of the cell. The choice of *δ* depends on the target function *f* and semi-metric *`*. More detail on that can be found in the next paragraph. Along which dimension to split is determined from the tree-depth ~ at which the cell is situated.

**Example 3.3.** Figure 3.6 shows the graph of two example functions we seek to maximize on the domain *X* = [0*,* 1]. Notably, the Garland function is not differentiable at some points on the domain and has no Lipschitz constant. The function is also not differentiable at the optimizer *π/*6. But there exists a semi-metric for which the Garland function is locally smooth around the optimizer.

**(a)** Sine and quadratic: *f*1(*x*) = 0*.*25 sin(50*x*) · sin(10*x*) − (*x* − 0*.*75)<sup>2</sup> with a scaled Euclidean metric fitted to the optimizer.

**(b)** Garland function: *f*2(*x*) = 4*x* (1 − *x*)( <sup>3</sup> <sup>4</sup> + 1 4 (1 − p |sin(60*x*)|)) and the semi-metric *`*(*x, y*) = *β*k*x* − *y*k <sup>1</sup>*/*<sup>2</sup> fitted to the optimizer.

**Figure 3.6:** Example functions that are locally smooth around the optimized for a semi-metric *`*.

³Symbols with a crossing bar as in ~ are used to denote height-indices in a tree.

To show convergence of DOO, the following assumptions are made for *f* and its domain [Mun+14].

For every region (~*, i*) containing the optimizer x <sup>∗</sup> ∈ *X*~*,i*, we have *f*(x~*,i*) + *δ*(~) ≥ *f*(x~*,i*) + *`*(x~*,i,* x ∗ ) ≥ *f* ∗ . Since the leaf nodes always cover the entire function domain, cells with x~*,i* + *δ*(~) *< f* <sup>∗</sup> are never expanded as they are dominated by the leaf node containing the optimizer. So the cells potentially expanded at depth ~ are *I*<sup>~</sup> = {(~*, i*) : *f*(x<sup>~</sup>*,i*) + *δ*(~) ≥ *f* <sup>∗</sup>}.

**(Stochastic) Simultaneous Optimistic Optimization** DOO requires no global Lipschitz bound of the target function. But it requires knowledge of a semi-metric *`* that is smooth around the optimizer. The *`* is not known in many cases. The Simultaneous Optimistic Optimization (SOO) algorithm [Mun11] adopts ideas from the DIRECT algorithm [JPS93] to achieve nearly the same convergence results as DOO even without knowledge of *`*. SOO still assumes the existence of a such a semi-metric *`*. But it suffices to show the existence of any such semi-metric for the convergence analysis without actually using it in the algorithm. In many cases, the function *f* itself can be used to construct a suitable semi-metric! Take any norm k · k for the function domain *X*. The semi-metric *`*(x*,* y) for points x*,* y on *X* is constructed as follows. With *η* = kx − yk the distance on the

domain norm, the distance on *`* is the difference between the optimizer and the worst point in the *η*-ball around the optimizer:

$$\tilde{\ell}(\eta) = \sup\_{\|x - \mathfrak{a}^\*\| \le \eta} (f^\* - f(\mathfrak{x})), \quad \ell(\mathfrak{x}, \mathfrak{y}) = \tilde{\ell}(\|x - \mathfrak{y}\|) \tag{3.8}$$

We forego full convergence proofs at this point and refer to the original paper [Mun11]. Stochastic SOO (StoSOO) [VCM13] extends SOO to the optimization of stochastic functions. Cells are split only after *κ* evaluations. The mean of the sampled values at the cell midpoint is used for the evaluation. The authors of [VCM13] provide a convergence bound for the expected regret on the order of *O*(log<sup>2</sup> (*ι*)*/* √ *ι*) where *ι* is the number of samples taken from the stochastic function *f*. The tuning parameter *η* controls how much emphasis the algorithm is putting on exploration, i.e. the tradeoff between exploration and exploitation.

Algorithm 11 shows a simplified version of StoSoo for the *n*-dimensional unit cube where cells are split into three child cells after *κ* evaluations. We will now explain the major changes compared to DOO. First, the visited nodes are split into leaf nodes *L* and internal nodes *N*. The internal nodes have been sampled *κ* times and are no longer evaluated. Second, since we do know the semi-metric *`* (and hence the cell diameter *δ*), an upper confidence bias is added to the cell evaluation for the selection. Third, the algorithm iterates over the depth-level ~ of the cells. One cell is selected at every level (if there is an improvement compared to the previous levels) and the cell is either sampled again or split.

### **Planning for Parameterized Action Sequences**

To integrate StoSOO with MCTS, we introduce so-called hybrid trees. Hybrid trees contain nodes for actions and parameters. Hybrid trees are bipartite as a parameter selection must follow an action selection and vice versa. Partiallyobservable hybrid trees (POHT) use three types of nodes: actions, parameters and the resulting observations. See Figure 3.7 for an example. Note that hybrid trees not only grow at the leaf nodes. The paramter-nodes represent a cell in the parameter space Θ*<sup>a</sup>* of the associated action *a*. As the cells of a continuous domain can be partitioned indefinitely often, a hybrid tree can grow new branches at the parameter-nodes during planning.

**Algorithm 11** Stochastic Simultaneous Optimistic Optimization (StoSOO) on the [0*,* 1]*<sup>n</sup>* cube. Cells are split after having been sampled *κ* times.

1: **procedure** StoSOO(*f, κ, ι, η*) 2: *L* ← {(0*,* 1)}*,* x[0*,* 1] ← 1 1 2 3: *q*[0*,* 1] ∼ *f*(x[0*,* 1])*, n*[0*,* 1] ← 1 4: **while** less than *ι* samples taken **do** 5: *q* max ← −∞ 6: **for** ~ = 0*, . . . ,* min(depth(*L*)*,* ~ max) **do** 7: *<sup>i</sup>* <sup>←</sup> arg max*<sup>j</sup>*:(~*,j*)∈*<sup>L</sup> <sup>q</sup>*[~*, j*] + <sup>q</sup>log(*<sup>ι</sup>* <sup>2</sup>*/η*) 2*n*[~*,j*] 8: *<sup>r</sup>* <sup>←</sup> *<sup>q</sup>*[~*, i*] + <sup>q</sup>log(*<sup>ι</sup>* <sup>2</sup>*/η*) 2*n*[~*,i*] 9: **if** *r* ≥ *q* max **then** 10: *q* max ← *r* 11: **if** *n*[~*, i*] *< κ* **then** 12: SampleStoSOO(~*, i*) 13: **else** 14: SplitStoSOO(~*, i*) 15: **return** arg max(~*,i*)*,n*[~*,i*]*>*<sup>0</sup> *q*[~*, i*] 1: **procedure** SampleStoSOO(~*, i*) 2: *y* ∼ *f*(x[~*, i*]) 3: *n*[~*, i*] ← *n*[~*, i*] + 1 4: *q*[~*, i*] ← *q*[~*, i*] + *<sup>y</sup>*−*q*[~*,i*] *n*[~*,i*] 1: **procedure** SplitStoSOO(~*, i*) 2: *v* ← max{*l* : x[~ + 1*, l*] 6= ∅} 3: *d* ← mod(~*, n*) + 1 4: o ← ν*<sup>d</sup>* 1 3 <sup>b</sup>~*/n*c+1 5: **for** *j* ∈ {1*, . . . ,* 3} **do** 6: *L* ← *L* ∪ {(~ + 1*, v* + *j*)} 7: x[~ + 1*, v* + *j*] ← x[~*, i*] + (*j* − 2)o 8: SampleStoSOO(~ + 1*, v* + *j*) 9: *L* ← *L* \ {(~*, i*)}

**Figure 3.7:** Partially Observable Hybrid Tree. The rectangle encompassing several circular nodes denotes a subtree for the optimisation of an action parameter with optimistic optimization. The tree is not fully explored for better visual representation.

For planning in POHT, we want to combine OO with MCTS. The upper confidence bound is used to select discrete actions and OO is used to select and iteratively refine the choice of action parameters. The StoSOO algorithm samples from the stochastic target function *f* several times within the SplitStoSOO procedure. This prevents its unmodified use for simulation-based planning – without backtracking – in sequential decision-making settings. Instead, we want every call to the OO subproblem to return exactly one parameter combination *θ* to continue the playout with the following steps in the scenario. The accumulated reward is then used to update the statistics for the involved branches of the scenario tree.

Previous authors have used OO for sequential decision-making with a continuous action-space [MWL11; Bus+13; BPM18]. The TrailBlazer algorithm of [GVM16] combines discrete action selection with continuous search. But it uses backtracking to return to a previous position in the scenario tree. In this thesis, we

want to avoid storing the system state for backtracking search. So the algorithms can just as well be performed with playouts in physical experiments only.

We now develop the novel Partially Observable Hybrid Tree Planning (POHTP) algorithm that combines MCTS – and in particular the approach for partiallyobservable planning from to the POMCP algorithm – with Optimistic Optimization for parameter selection on continuous domains. See Algorithm 12 for the full details. In every step for sequential decision making, the algorithm is presented with the choice of several discrete actions that are each parameterized from a continuous domain. If there is only one action, POHTP reduces to StoSOO as a special case. On the other hand, if the actions have no parameters, POHTP reduces to a variation of POMCP. The difference to the original POMCP is the update mechanism that maximizes over possible choices to compute the V-value estimate (instead of taking the empirical reward from the previous plays), full playouts and a switch between exploitation and exploration within each playout that is explained in the next paragraph.