# Optimization with Temporal and Logical Specifications
via Generalized Mean-based Smooth Robustness Measures

Samet Uzun<EMAIL_ADDRESS>Purnanand Elango<EMAIL_ADDRESS>Pierre-Loïc Garoche
<EMAIL_ADDRESS>Behçet Açıkmeşe<EMAIL_ADDRESS>∗William E. Boeing
Department of Aeronautics & Astronautics, University of Washington, Seattle,
WA, 98195  †École Nationale de l’Aviation Civile, Université de Toulouse,
France

###### Abstract

This paper introduces a generalized mean-based $\mathcal{C}^{1}$-smooth
robustness measure over discrete-time signals (D-GMSR) for signal temporal
logic (STL) specifications. In conjunction with its
$\mathcal{C}^{1}$-smoothness, D-GMSR is proven to be both sound and complete.
Furthermore, it demonstrates favorable gradient properties and addresses
locality and masking problems, which are critical for numerical optimization.
The $\mathcal{C}^{1}$-smoothness of the proposed formulations enables the
implementation of robust and efficient numerical optimization algorithms to
solve problems with STL specifications while preserving their theoretical
guarantees. The practical utility of the proposed robustness measure is
demonstrated on two real-world trajectory optimization problems: i) quadrotor
flight, and ii) autonomous rocket landing. A sequential convex programming
(SCP) framework, incorporating a convergence-guaranteed optimization algorithm
(the prox-linear method) is used to solve inherently non-convex trajectory
optimization problems with STL specifications. The implementation is available
at https://github.com/UW-ACL/D-GMSR

###### keywords:

Signal Temporal Logic; Robustness Measure; Sequential Convex Programming;
Trajectory Optimization; Optimal Control;

captionUnknown document class ††thanks: This research was supported by AFOSR
grant FA9550-20-1-0053 and ONR grant N00014-20-1-2288; Government sponsorship
is acknowledged.
Corresponding author: Samet Uzun. ∗, ∗, †, ∗

##  1 Introduction

Temporal logic specifications [1] provide a systematic approach for encoding
high-level specifications and complex system behaviors in mission planning and
control. The application of temporal logic specifications in dynamical systems
can be found in autonomous vehicles [2], high-dimensional manipulators [3],
and multi-agent systems [4], including mobile robots [5, 6], and UAVs [7, 8].
Temporal logic specifications are also integrated into learning algorithms to
ensure that results satisfy certain properties in prediction [9] and policy
synthesis [10, 11]. Notably, useful formulations have emerged to facilitate
systematic incorporation of such constraints [12].

Linear temporal logic (LTL) is one of the most popular temporal logic
frameworks for expressing design requirements for discrete event systems [13].
Finite state automata [14] have been traditionally employed to control
discrete-time dynamical systems with LTL specifications. A receding horizon
framework is proposed in [15] to alleviate the computational complexity of
synthesizing finite state automata, along with the associated software toolbox
[16]. [17] introduced a general technique, formulating LTL specifications as
mixed-integer linear constraints, enabling optimization algorithms to address
these problems. Metric temporal logic [18] and metric interval temporal logic
[19] serve as natural extensions of LTL, to extend the discrete-time
specifications to continuous time. Signal temporal logic (STL) combines the
dense temporal modalities of metric interval temporal logic (MITL) with
predicates specifying properties over real-valued signals, motivated by the
application of verifying and monitoring temporal properties of continuous-time
signals in continuous and hybrid systems [20]. A robust semantic for MITL,
which generates a real number to quantify the degree of specification
satisfaction instead of a Boolean evaluation, is introduced in [21]. A similar
robustness measure, called space robustness (SR), and a time robustness
concept aimed at improving robustness against timing uncertainties within STL
is introduced in [22]. SR utilizes $\min$ and $\max$ functions,  which we
refer to as key functions of SR, for logical and temporal operators of STL.
Hence, the introduced robustness measures allow the encoding of temporal
specifications as cost functions, transforming the synthesis problem into an
optimization problem [23, 24]. While SR is initially formulated for
continuous-time signals [22], a discrete-time counterpart (D-SR), which
resembles the robustness measure proposed for MITL in [21], has been explored
in related literature [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. The
inherent design of the introduced robustness measures leads to solutions in
the form of mixed-integer programming (MIP) problems for both D-SR [25, 26]
and time robustness [27, 28]. The formulation of D-SR with the $\min$ and
$\max$ functions has been addressed by smoothing these functions using the
well-known log-sum-exponential (LSE) functions in [29]. The work is further
extended in [30] to address the notion of cumulative robustness. The LSE
approximation is asymptotically sound and asymptotically complete, serving as
an under-approximation for the $\min$ function and an over-approximation for
the $\max$ function. The LSE approximation for the $\max$ function is replaced
with an under-approximation in [31], resulting in a sound, asymptotically
complete and smooth space robustness (D-SSR) for STL. A polynomial smooth
approximation is used in [32] for a more accurate approximation of the $\min$
and $\max$ functions. Although the number of variables of these functions that
approximates the $\min$ and $\max$ functions cannot be greater than two,
recursive modeling of STL constraints in this paper causes all $\min$ and
$\max$ functions to perform on only two variables. While the robust semantics
introduced for STL specifications demonstrate utility in verification and
monitoring [22], the smooth versions may pose challenges for optimization,
primarily related to concerns about locality and masking issues defined in
[36]. Locality is defined as the dependency of the robustness measure on
signals at a single time instant, whereas masking describes the situation
where the robustness measure is solely determined by one of the subformulas
combined with conjunction and disjunction operators. Utilizing an averaging
approach that considers the entire signal and subformulas, instead of relying
on specific time instances or subformulas through the $\min$ and $\max$
functions, proves beneficial in mitigating these issues. Averaged-STL is
introduced in [37] to address the falsification problem in continuous-time
signals. This approach involves creating a weighted time average over specific
STL operators. However, it does not eliminate the use of the $\min$ and $\max$
functions, resulting in the continuation of non-smooth STL semantics. In [33],
discrete average space robustness (DASR) and its simplified version (DSASR)
are introduced as computationally more tractable alternatives to D-SR.
However, both DASR and DSASR incorporate $\min$ and $\max$ functions, leading
to non-smoothness in certain STL operators. Additionally, ensuring the
soundness of the STL formula necessitates imposing additional constraints. The
arithmetic-geometric mean robustness (AGM), presented in [34], emphasizes both
the frequency of satisfaction of STL specifications and the robustness of each
specification at every time instant. The AGM also addresses the locality and
masking problems. However, the introduced robustness measure switches between
arithmetic and geometric means with if-else conditions to satisfy the STL
specification, leading to non-smoothness. Additionally, the robustness degree
can only cover a specific set of STL specifications. Both D-SR and AGM are
generalized in [35] by introducing weights associated with STL operators to
enable the definition of user preferences for priorities and importance.

This paper introduces a generalized mean-based $\mathcal{C}^{1}$-smooth
robustness measure over discrete-time signals (D-GMSR). D-GMSR is constructed
by replacing the key functions of D-SR, $\min$, and $\max$ functions, with two
$\mathcal{C}^{1}$-smooth functions that are specific compositions of
generalized means. These functions have the same sign with $\min$ and $\max$
functions, respectively. Therefore, D-GMSR is both sound and complete, making
it reliable and non-conservative. Furthermore, the first-order information of
these functions satisfies some desirable properties critical for numerical
optimization purposes and addresses the issues of locality and masking. Table
1 compares D-GMSR with existing robustness measures from the literature. The
$\mathcal{C}^{1}$-smoothness of the robustness measure facilitates the
application of a broad range of numerical solvers for optimization problems
with STL specifications while preserving their theoretical guarantees. We
provide some illustrative numerical examples of D-GMSR and demonstrate its
practical utility in trajectory optimization for quadrotor flight and
autonomous rocket landing problems solved via sequential convex programming
(SCP) framework [38, 39, 40, 41]. The implementation is available at
https://github.com/UW-ACL/D-GMSR

| [22] | [29] | [31] | [34] | D-GMSR
---|---|---|---|---|---
Smoothness | $\times$ | ✓ | ✓ | $\times$ | ✓
Soundness | ✓ | $\circ$ | ✓ | ✓ | ✓
Completeness | ✓ | $\circ$ | $\circ$ | ✓ | ✓
Locality & Masking | $\times$ | $\times$ | $\times$ | ✓ | ✓

Table 1: Comparison with the previous robustness measures
($\circ$: Satisfied Asymptotically)

The paper is organized as follows. Section 2 introduces the syntax of STL and
smooth space robustness, quantifying STL specification satisfaction, yet
highlighting its limitations resulting in undesirable optimization outcomes
where the specification is not met. These limitations are addressed by the
proposed generalized mean-based $\mathcal{C}^{1}$-smooth robustness measure
(D-GMSR) in Section 3. Section 4 demonstrates the practical utility of D-GMSR
in trajectory optimization for quadrotor flight and rocket landing problems
solved within an SCP framework. The paper is concluded in Section 5.

Notation: The real numbers are denoted as $\mathbb{R}$, where $\mathbb{R}_{+}$
and $\mathbb{R}_{++}$ represent the sets of non-negative and positive real
numbers, respectively. Similarly, $\mathbb{Z}_{+}$ and $\mathbb{Z}_{++}$
represent the sets of non-negative and positive integers. $\mathbb{R}^{n\times
m}$ and $\mathbb{R}^{n}$ denote the sets of real $n\times m$ matrices and
$n\times 1$ vectors, respectively. $\mathrm{diag}(A)\in\mathbb{R}^{n}$
represents the diagonal elements of a matrix $A\in\mathbb{R}^{n\times n}$.
$\bm{0}$ denotes a zero vector, while $\bm{1}$ represents a vector of ones in
appropriate dimensions. The vector $e_{i}$ is defined such that its $i$th
element is $1$ and all other elements are $0$. The cardinality of set $S$ is
denoted by $\mathrm{card}(S)$. $\top$ and $\bot$ represent boolean true and
false, respectively. The space of continuous and $k$-times differentiable
functions on a domain $X$ is denoted by $\mathcal{C}^{k}(X)$, with members
referred to as $\mathcal{C}^{k}$-smooth. A vector $x\in\mathbb{R}^{n}$ is
represented as $(x_{1},x_{2},\dots,x_{n})$.
$[a{:}b]:=\\{k\in\mathbb{Z}_{+}:a\leq k\leq b;a,b\in\mathbb{Z}_{+}\\}$.
$x_{[a{:}b]}:=(x_{a},x_{a+1},\dots,x_{b})$. The functions $\min(x)$ and
$\max(x)$ return the minimum and maximum elements of the $x\in\mathbb{R}^{n}$
vector, respectively.
$|x|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}$ is
defined as either $|x|_{+}$ or $|x|_{-}$.
$|x|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{p}=(f((x_{1},0))^{p},\dots,f((x_{n},0))^{p})$,
where $f=\max$ if
$\text{\framebox{\rule{0.0pt}{4.30554pt}\rule{4.30554pt}{0.0pt}}}=+$, and
$\min$ otherwise.

##  2 Preliminaries

In this section, we provide a concise overview of the syntax and a robust
semantic of STL tailored for discrete-time signals.

###  2.1 Signal Temporal Logic (STL)

Let $x:\mathbb{Z}_{+}\to\mathbb{R}^{n}$ be a discrete-time signal and $\mu$ be
a Boolean-valued predicate such that

$\displaystyle\mu:=$ $\displaystyle\;(f(x)\geq 0)$ $\displaystyle=$
$\displaystyle\begin{cases}\top&\text{if $f(x)\geq 0$}\\\
\bot&\text{otherwise}.\end{cases}$

where $f:\mathbb{R}^{n}\to\mathbb{R}$ is a real-valued function referred to as
the predicate function. The following assumption is needed for numerical
optimization purposes.

######  Assumption 1.

All predicate functions are $\mathcal{C}^{1}$-smooth.

Let $I$ be a non-empty time interval such that $I=[a{:}b]$. The syntax of STL
is defined recursively as follows:

$\varphi::=\mu\;|\;\neg\varphi_{1}\;|\;\varphi_{1}\wedge\varphi_{2}\;|\;\varphi_{1}\bm{U}_{I}\varphi_{2}$
(1)

where $\varphi$ is the STL specification (or STL formula), $\varphi_{1}$ and
$\varphi_{2}$ are the STL sub-specifications (or STL sub-formulas) , $\neg$
and $\wedge$ are the Boolean negation and conjunction, respectively, and
$\bm{U}_{I}$ is the until temporal operator over bounded interval $I$. In
order to construct an STL specification, different predicates or STL sub-
formulas are recursively combined using the Boolean negation, conjunction, and
until operators as in the examples in Table 2.

STL specification | Formula
---|---
disjunction | $\varphi_{1}\vee\varphi_{2}$ | $\neg(\neg\varphi_{1}\wedge\neg\varphi_{2})$
implication | $\varphi_{1}\implies\varphi_{1}$ | $\neg\varphi_{1}\vee\varphi_{2}$
eventually | $\bm{F}_{I}\varphi$ | $\top\bm{U}_{I}\varphi$
always | $\bm{G}_{I}\varphi$ | $\neg\bm{F}_{I}\neg\varphi$
Table 2: Construction of some STL specifications

Formally, the semantics of an STL formula $\varphi$ defines what it means for
a signal $x$ to satisfy $\varphi$ at time point $k$, denoted as
$(x,k)\vDash\varphi$. Conversely, if the specification is not satisfied, it is
denoted as $(x,k)\nvDash\varphi$.

######  Remark 1.

The $\bm{U}_{I}$ operator in the STL syntax (1), and consequently the
$\bm{F}_{I}$ and $\bm{G}_{I}$ operators, can be represented over discrete-time
signals by utilizing the first three elements, namely the predicate, negation,
and conjunction.

###  2.2 Space Robustness for STL

While determining if the signal meets the STL specification with a Boolean
evaluation, the robust semantics yield a real number that quantifies the
degree of satisfaction of a specification. A space robustness (SR) concept, a
robustness measure for STL, is presented in [22] by utilizing $\min$ and
$\max$ functions for logical and temporal operators of STL. In this approach,
temporal specifications can be represented as cost functions, transforming a
synthesis problem into an optimization problem [23, 24]. While SR is initially
developed for continuous-time signals, its discrete-time version (D-SR) has
been explored in various studies [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35].

######  Definition 2.1.

(D-SR [22, 25]) The space robustness [22] defines a real-valued function
$\rho^{\varphi}$ of signal $x$ and time $k$ such that
$(x,k)\models\varphi\iff\rho^{\varphi}(x,k)\geq 0$. The complete D-SR is
defined as follows:

$\displaystyle\rho^{\mu}(x,k):=f(x_{k})$
$\displaystyle\rho^{\neg\varphi}(x,k):=-\rho^{\varphi}(x,k)$
$\displaystyle\rho^{\varphi_{1}\wedge\varphi_{2}}(x,k):=\min((\rho^{\varphi_{1}}(x,k),\rho^{\varphi_{2}}(x,k)))$
$\displaystyle\rho^{\varphi_{1}\vee\varphi_{2}}(x,k):=\max((\rho^{\varphi_{1}}(x,k),\rho^{\varphi_{2}}(x,k)))$
$\displaystyle\rho^{\varphi_{1}\implies\varphi_{2}}(x,k):=\max((-\rho^{\varphi_{1}}(x,k),\rho^{\varphi_{2}}(x,k)))$
$\displaystyle\rho^{\bm{F}_{[a:b]}\varphi}(x,k):=\max((\rho^{\varphi}(x,k+a),\rho^{\varphi}(x,k+a+1),$
$\displaystyle\hskip 93.89418pt\dots,\rho^{\varphi}(x,k+b)))$
$\displaystyle\rho^{\bm{G}_{[a:b]}\varphi}(x,k):=\min((\rho^{\varphi}(x,k+a),\rho^{\varphi}(x,k+a+1),$
$\displaystyle\hskip 93.89418pt\dots,\rho^{\varphi}(x,k+b)))$
$\displaystyle\rho^{\varphi_{1}\bm{U}_{[a:b]}\varphi_{2}}(x,k):=\max((\min((\rho^{\bm{G}_{[a:a]}\;\varphi_{1}}(x,k),$
$\displaystyle\hskip 135.43495pt\rho^{\varphi_{2}}(x,k+a))),$
$\displaystyle\hskip
109.25847pt\min((\rho^{\bm{G}_{[a:a+1]}\;\varphi_{1}}(x,k),$
$\displaystyle\hskip 135.43495pt\rho^{\varphi_{2}}(x,k+a+1))),$
$\displaystyle\hskip 118.79033pt\vdots$ $\displaystyle\hskip
109.25847pt\min((\rho^{\bm{G}_{[a:b]}\;\varphi_{1}}(x,k),$
$\displaystyle\hskip 135.43495pt\rho^{\varphi_{2}}(x,k+b)))))$

The $\min$ and $\max$ functions are key functions of D-SR, and any STL
specification can be formulated by composing the predicate functions with
these key functions. However, the non-smoothness of the $\min$ and $\max$
functions results in optimization problems involving STL specifications
becoming MIP problems [25, 26]. In [31], a smooth version of D-SR (D-SSR) is
proposed by replacing $\min$ and $\max$ functions with their under-
approximations as follows:

$\displaystyle\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(x)$
$\displaystyle:=-\frac{1}{\kappa}\mathrm{log}\bigg{(}\sum_{i=1}^{n}e^{-\kappa
x_{i}}\bigg{)}$ (2a)
$\displaystyle\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(x)$
$\displaystyle:=\frac{\sum_{i=1}^{n}x_{i}e^{\kappa
x_{i}}}{\sum_{i=1}^{n}e^{\kappa x_{i}}}$ (2b)

where $x\in\mathbb{R}^{n}$ and $\kappa>0$ is a sharpness parameter determining
the accuracy of the approximation. Thus, a smooth, sound, and asymptotically
complete robustness measure $\tilde{\rho}^{\varphi}(x,k)$ is achieved for STL,
enabling the utilization of gradient-based algorithms to solve optimization
problems involving STL specifications.

###  2.3 Gradient properties of the D-SSR

In D-SR, the values of multiple predicate functions at various time instances
can be composed using the $\min$ and $\max$ functions, such as

$\max((\min((f_{1}(x_{1}),f_{2}(x_{1}))),f_{1}(x_{2}))).$

Due to the nature of $\min$ and $\max$ functions, the robustness value of an
STL specification is typically determined by the value of a specific predicate
function at a single time instant. The values of all other predicate functions
and values of predicate functions at other time instances are ignored, which
are called locality and masking, respectively, as introduced in [36].
Consequently, the derivatives of
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
and
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
functions with respect to one of the input variables are nearly zero if the
minimum or maximum value are not attained by this variable among others.

To demonstrate this, suppose that $f_{i}:\mathbb{R}^{n}\to\mathbb{R}$ is a
function of $x\in\mathbb{R}^{n}$ for all $i\in[1{:}m]$,
$S_{\textrm{min}}=\\{i:f_{i}(x)\leq f_{j}(x),\forall j\in[1{:}m]\\}$, and
$S_{\textrm{max}}=\\{i:f_{i}(x)\geq f_{j}(x),\forall j\in[1{:}m]\\}$. For
sufficiently large $\kappa$, if $y=(f_{1}(x),f_{2}(x),\dots,f_{m}(x))$, then

$\displaystyle\nabla_{x}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(y)$
$\displaystyle=\sum_{i=1}^{m}a_{i}\nabla_{x}f_{i}(x)\text{ and}$
$\displaystyle\nabla_{x}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(y)$
$\displaystyle=\sum_{i=1}^{m}b_{i}\nabla_{x}f_{i}(x),\text{ where}$
$\displaystyle a_{i}$
$\displaystyle\approx\begin{cases}\frac{1}{\mathrm{card}(S_{\textrm{min}})}&\text{if
}i\in S_{\textrm{min}}\\\ 0&\text{otherwise}\end{cases},$ $\displaystyle
b_{i}$
$\displaystyle\approx\begin{cases}\frac{1}{\mathrm{card}(S_{\textrm{max}})}&\text{if
}i\in S_{\textrm{max}}\\\ 0&\text{otherwise}\end{cases}.$

Therefore, using D-SSR to account for an STL specification within a gradient-
based iterative optimization algorithm leads to the consideration of only the
gradient of one predicate function at a single time instant, $f_{i}(x_{k})$.

######  Example 2.2.

In Figure 1, we illustrate the variations in the predicate function values
during the maximization of
$\tilde{\rho}^{\varphi_{1}\wedge\varphi_{2}\wedge\dots\wedge\varphi_{5}}(x,k)$
and $\tilde{\rho}^{\varphi_{1}\vee\varphi_{2}\vee\dots\vee\varphi_{5}}(x,k)$
functions, respectively. The sharpness parameter determining the accuracy of
the approximation of the $\min$ and $\max$ functions (2) is set to
$\kappa=25$. Here, $\tilde{\rho}^{\varphi_{i}}(x,k)=f_{i}(x_{k})\in\mathbb{R}$
and $f_{i}:\mathbb{R}\to\mathbb{R}$ for all $i\in[1{:}5]$. Therefore,

$\displaystyle\tilde{\rho}^{\varphi_{1}\wedge\varphi_{2}\wedge\dots\wedge\varphi_{5}}(x,k)$
$\displaystyle=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(y),$
$\displaystyle\tilde{\rho}^{\varphi_{1}\vee\varphi_{2}\vee\dots\vee\varphi_{5}}(x,k)$
$\displaystyle=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(y),$

where $y=(f_{1}(x_{k}),f_{2}(x_{k}),\dots,f_{5}(x_{k}))$. Alternatively,
maximizing $\tilde{\rho}^{\bm{G}_{[1:5]}\varphi}(x,0)$ and
$\tilde{\rho}^{\bm{F}_{[1:5]}\varphi}(x,0)$ functions, respectively, where
$\tilde{\rho}^{\varphi}(x,k)=f(x_{k})\in\mathbb{R}$ for all $k\in[1{:}5]$ and
$f:\mathbb{R}\to\mathbb{R}$, leads to similar changes in predicate function
values across different time instances since

$\displaystyle\tilde{\rho}^{\bm{G}_{[1:5]}\varphi}(x,0)$
$\displaystyle=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(y),$
$\displaystyle\tilde{\rho}^{\bm{F}_{[1:5]}\varphi}(x,0)$
$\displaystyle=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(y),$

where $y=(f(x_{1}),f(x_{2}),\dots,f(x_{5}))$.

During the maximization of
$\tilde{\rho}^{\varphi_{1}\wedge\varphi_{2}\wedge\dots\wedge\varphi_{5}}(x,k)$
function, only the value of the predicate function that attains the minimum
value increases notably with each iteration and the function with the minimum
value changes after a certain number of iterations. Hence, the values of all
predicate functions become positive after a large number of iterations. On the
other hand, during the maximization of
$\tilde{\rho}^{\varphi_{1}\vee\varphi_{2}\vee\dots\vee\varphi_{5}}(x,k)$
function, only the value of the predicate function that attains the maximum
value increases notably with each iteration. Therefore, the predicate function
with the highest value does not change and the relevant function becomes
positive after a large number of iterations. These characteristics of
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
and
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
functions may lead to undesirable results, i.e., the STL specification is not
met, as illustrated by an example in Section 4.

These undesirable characteristics might be alleviated by introducing a small
smoothing parameter $\kappa$ for the $\min$ and $\max$ functions in the
initial iterations of the numerical optimization algorithm, gradually
increasing $\kappa$ with subsequent iterations to achieve a more accurate
approximation similar to the homotopy approach employed in the SCP algorithm
[42]. However, this approach requires the inclusion of additional decision
variables in the algorithm and may potentially disrupt its useful convergence
guarantees.

Figure 1:  The changes in the values of $f_{i}(x_{k})$ (or $f(x_{i})$) and the
corresponding robustness values while maximizing the conjunction (or always)
and disjunction (or eventually) operators using D-SSR are presented on the
left and right sides of the figure, respectively. The color gradient, ranging
from purple to red in the figures above, represents the values of the
functions from $f_{1}(x_{k})$ to $f_{5}(x_{k})$ (or $f(x_{1})$ to $f(x_{5})$),
respectively.

##  3  Generalized Mean-based Smooth Robustness over Discrete-time Signals
(D-GMSR)

To remedy the issues in previously proposed robustness measures, this section
introduces generalized mean-based $\mathcal{C}^{1}$-smooth robustness measure
over discrete-time signals (D-GMSR). We first define functions that are the
key for D-GMSR and then introduce D-GMSR by replacing $\min$ and $\max$
functions of D-SR in Definition 2.1 with the proposed key functions. Here we
also provide proofs of the $\mathcal{C}^{1}$-smoothness, soundness, and
completeness properties of the D-GMSR. We then present an analysis of the
gradients of the proposed key functions and provide some numerical
illustrative examples for the application of D-GMSR.

We propose the following $\mathcal{C}^{1}$-smooth function, constructed
through generalized means, characterized by the same sign as the $\min$
function. This proposed function will replace the $\min$ function of D-SR in
Definition 2.1.

${}^{\wedge}h_{p,w}^{\epsilon}(y):=\bigg{(}M_{0,w}^{\epsilon}(|y|_{+}^{2})\bigg{)}^{\frac{1}{2}}-\bigg{(}M_{p,w}^{\epsilon}(|y|_{-}^{2})\bigg{)}^{\frac{1}{2}},$
(3)

where

$\displaystyle M_{0,w}^{\epsilon}(z)$
$\displaystyle:=\Big{(}\epsilon^{\bm{1}^{T}w}+\prod_{i=1}^{n}z_{i}^{w_{i}}\Big{)}^{1/\bm{1}^{T}w}$
(4a) $\displaystyle M_{p,w}^{\epsilon}(z)$
$\displaystyle:=\Big{(}\epsilon^{p}+\frac{1}{\bm{1}^{T}w}\sum_{i=1}^{n}w_{i}z_{i}^{p}\Big{)}^{1/p}$
(4b)

$y\in\mathbb{R}^{n}$, $z\in\mathbb{R}^{n}_{+}$, $\epsilon\in\mathbb{R}_{++}$,
$p\in\mathbb{Z}_{++}$, $w\in\mathbb{Z}_{++}^{n}$. We define another key
function that will replace the $\max$ function of D-SR in Definition 2.1 as
${}^{\vee}h_{p,w}^{\epsilon}(y):=-{}^{\wedge}h_{p,w}^{\epsilon}(-y)$.

######  Remark 2.

If $\epsilon=0$, then $M_{0,w}^{\epsilon}(z)$ and $M_{p,w}^{\epsilon}(z)$ are
weighted geometric mean and weighted power mean of $z$, respectively.
Furthermore,
$M_{0,w}^{\epsilon},M_{p,w}^{\epsilon}:\mathbb{R}^{n}_{+}\to[\epsilon,\infty)$,
and
$M_{0,w}^{\epsilon},M_{p,w}^{\epsilon}\in\mathcal{C}^{\infty}(\mathbb{R}_{+}^{n})$.
The values of the ${}^{\wedge}h_{p,w}^{\epsilon}$ and
${}^{\vee}h_{p,w}^{\epsilon}$ functions vary depending on the non-negativity
and non-positivity, respectively, of their input vectors as follows:
$\begin{aligned}
{}^{\wedge}h_{p,w}^{\epsilon}(y)&=\begin{cases}(M_{0,w}^{\epsilon}(|y|_{+}^{2}))^{0.5}-\epsilon^{0.5}&\text{if
$y_{i}\geq 0$, $\forall i\in[1{:}n]$ }\\\
\epsilon^{0.5}-(M_{p,w}^{\epsilon}(|y|_{-}^{2}))^{0.5}&\text{otherwise}\end{cases}\\\
{}^{\vee}h_{p,w}^{\epsilon}(y)&=\begin{cases}\epsilon^{0.5}-(M_{0,w}^{\epsilon}(|y|_{-}^{2}))^{0.5}&\text{if
$y_{i}\leq 0$, $\forall i\in[1{:}n]$ }\\\
(M_{p,w}^{\epsilon}(|y|_{+}^{2}))^{0.5}-\epsilon^{0.5}&\text{otherwise}.\end{cases}\end{aligned}$
The value of $\epsilon$, and the parameters associated with the generalized
means, $p$, and $w$, can be chosen by considering the underlying problem as
follows:

  * •

For any $\epsilon>0$, ${}^{\wedge}h_{p,w}^{\epsilon}$ and
${}^{\vee}h_{p,w}^{\epsilon}$ functions demonstrate
$\mathcal{C}^{1}$-smoothness, while preserving the soundness and completeness
properties of the D-GMSR (see Lemma 3.5). It is desirable to select an
arbitrarily small positive value for $\epsilon$ to ensure that the values of
$M_{0,w}^{\epsilon}$ and $M_{p,w}^{\epsilon}$ functions closely approximate
the weighted geometric mean and weighted power mean of their input.

  * •

While $M_{p,w}^{0}(z)$ represents a weighted arithmetic mean for $p=1$,
$(M_{p,w}^{\epsilon}(|z|_{+}^{2}))^{0.5}\to\max(z)$ and
$(M_{p,w}^{\epsilon}(|z|_{-}^{2}))^{0.5}\to\min(z)$ as $p\to\infty$. Hence,
users have the option to make the values of the
${}^{\wedge}h_{p,w}^{\epsilon}$ and ${}^{\vee}h_{p,w}^{\epsilon}$ functions as
similar to the values of the $\min$ and $\max$ functions under some conditions
(see Lemma 3.6).

  * •

The weight vector $w\in\mathbb{Z}_{++}^{n}$ allows users to express
preferences and priorities, indicating the relative importance of elements in
the vector $z\in\mathbb{R}^{n}_{+}$.

######  Definition 3.3.

(D-GMSR) D-GMSR is defined by replacing $\min$ and $\max$ functions of D-SR in
Definition 2.1 with ${}^{\vee}h_{p,w}^{\epsilon}(y)$ and
${}^{\wedge}h_{p,w}^{\epsilon}(y)$ functions, respectively. The complete
D-GMSR is defined as follows:  $\begin{aligned}
&\Gamma^{\mu}(x,k):=f(x_{k})\\\
&\Gamma^{\neg\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k):=-\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi}(x,k)\\\
&\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{1}\wedge\varphi_{2}}(x,k):={}^{\wedge}h_{p_{1},w_{1}}^{\epsilon_{1}}((\Gamma^{\varphi_{1}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k),\Gamma^{\varphi_{2}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)))\\\
&\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{1}\vee\varphi_{2}}(x,k):={}^{\vee}h_{p_{1},w_{1}}^{\epsilon_{1}}((\Gamma^{\varphi_{1}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k),\Gamma^{\varphi_{2}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)))\\\
&\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{1}\implies\varphi_{2}}(x,k):={}^{\vee}h_{p_{1},w_{1}}^{\epsilon_{1}}((-\Gamma^{\varphi_{1}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k),\Gamma^{\varphi_{2}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)))\\\
&\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\bm{F}_{[a{:}b]}\varphi}(x,k):={}^{\vee}h_{p_{1},w_{1}}^{\epsilon_{1}}((\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k+a),\\\
&\hskip 112.61578pt\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k+a+1),\\\
&\hskip
110.96556pt\dots,\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k+b)))\\\
&\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\bm{G}_{[a{:}b]}\varphi}(x,k):={}^{\wedge}h_{p_{1},w_{1}}^{\epsilon_{1}}((\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k+a),\\\
&\hskip 112.61578pt\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k+a+1),\\\
&\hskip
110.96556pt\dots,\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k+b)))\\\
&\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{1}\bm{U}_{[a{:}b]}\varphi_{2}}(x,k):={}^{\vee}h_{p_{1},w_{1}}^{\epsilon_{1}}(({}^{\wedge}h_{p_{2},w_{2}}^{\epsilon_{2}}((\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\bm{G}_{[a{:}a]}\;\varphi_{1}}(x,k),\\\
&\hskip
168.29785pt\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{2}}(x,k+a))),\\\
&\hskip
126.89948pt{}^{\wedge}h_{p_{2},w_{2}}^{\epsilon_{2}}((\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\bm{G}_{[a{:}a+1]}\;\varphi_{1}}(x,k),\\\
&\hskip
168.29785pt\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{2}}(x,k+a+1))),\\\
&\hskip 135.15059pt\vdots\\\ &\hskip
126.89948pt{}^{\wedge}h_{p_{2},w_{2}}^{\epsilon_{2}}((\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\bm{G}_{[a{:}b]}\;\varphi_{1}}(x,k),\\\
&\hskip
168.44026pt\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{2}}(x,k+b)))))\end{aligned}$
where $\bm{\epsilon}$, $\bm{p}$, and $\bm{w}$ represent sets consisting of an
arbitrary number of $\epsilon_{i}\in\mathbb{R}_{++}$,
$p_{i}\in\mathbb{Z}_{++}$, and $w_{i}\in\mathbb{Z}_{++}^{n}$, respectively,
required for the construction of ${}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}}$
and ${}^{\vee}h_{p_{i},w_{i}}^{\epsilon_{i}}$ functions, which in turn
construct the corresponding specification. For instance, let
$\varphi:=\bm{F}_{[1{:}5]}((f_{1}(x)\geq 0)\wedge(f_{2}(x)\geq 0))$. Then, the
cardinality of the sets $\bm{\epsilon}$, $\bm{p}$ and $\bm{w}$ used in
$\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi}(x,0)$ is $2$, where
$w_{1}\in\mathbb{Z}_{++}$ and $w_{2}\in\mathbb{Z}_{++}^{5}$.

###  3.1 Properties of the D-GMSR

In Lemma 3.4 and 3.5, we assume $\epsilon_{i}\in\mathbb{R}_{++}$,
$p_{i}\in\mathbb{Z}_{++}$, $w_{i}\in\mathbb{Z}_{++}^{n_{i}}$ for all
$i\in\mathbb{Z}_{++}$ and do not restate this for brevity.

######  Lemma 3.4.

($\mathcal{C}^{1}$-smoothness)
$\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(\cdot,k)\in\mathcal{C}^{1}(\mathbb{R}^{n})$
for any finite time STL specification $\varphi$ and time step
$k\in\mathbb{Z}_{+}$.

Proof: Given that Assumption 1 requires all predicate functions to be
$\mathcal{C}^{1}$-smooth, finite time STL specification $\varphi$ and time
step $k\in\mathbb{Z}_{+}$ if and only if each
${}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}}$ and
${}^{\vee}h_{p_{i},w_{i}}^{\epsilon_{i}}$ functions that construct the STL
specification $\varphi$ are $\mathcal{C}^{1}$-smooth on $\mathbb{R}^{n}$.
$|y|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{2}\in\mathbb{R}_{+}^{n}$
implies that
$M_{p_{i},w_{i}}^{\epsilon_{i}}(|y|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{2}),M_{0,w_{i}}^{\epsilon_{i}}(|y|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{2})\in[\epsilon,\infty]$.
Note that
$\nabla_{x}|x|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{2}=2|x|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}$
for $x\in\mathbb{R}^{n}$, which implies
$|\cdot|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{2}\in\mathcal{C}^{1}(\mathbb{R}^{n})$
Then
$M_{p_{i},w_{i}}^{\epsilon_{i}},M_{0,w_{i}}^{\epsilon_{i}}\in\mathcal{C}^{\infty}(\mathbb{R}_{+}^{n})$
implies that
$M_{p_{i},w_{i}}^{\epsilon_{i}}(|\cdot|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{2}),M_{0,w_{i}}^{\epsilon_{i}}(|\cdot|_{\text{\framebox{\rule{0.0pt}{3.01389pt}\rule{3.01389pt}{0.0pt}}}}^{2})\in\mathcal{C}^{1}(\mathbb{R}^{n})$.
Therefore,
${}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}},{}^{\vee}h_{p_{i},w_{i}}^{\epsilon_{i}}$
$\in\mathcal{C}^{1}(\mathbb{R}^{n})$.

######  Lemma 3.5.

(Soundness and Completeness) For any discrete-time signal $x$, time step
$k\in\mathbb{Z}_{+}$, and finite time STL specification $\varphi$,
$(x,k)\models\varphi\iff\Gamma^{\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)\geq
0$.

Proof: Given that D-SR in Definition 2.1 is both sound and complete [31], and
${}^{\vee}h_{p_{i},w_{i}}^{\epsilon_{i}}(y)=-{}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}}(-y)$,
analogous to $\max(y)=-\min(-y)$, proofs of the following assertions establish
the soundness and completeness of D-GMSR.

  * •

${}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}}(y)<0\iff\min(y)<0$,

  * •

${}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}}(y)>0\iff\min(y)>0$.

One can verify from (4) that

  * •

$M_{0,w_{i}}^{\epsilon_{i}}(|y|_{+}^{2})>\epsilon\iff y_{j}>0$ for all
$j\in[1{:}n]$,

  * •

$M_{p_{i},w_{i}}^{\epsilon_{i}}(|y|_{-}^{2})>\epsilon\iff\exists j\in[1{:}n]$
such that $y_{j}<0$,

  * •

$M_{0,w_{i}}^{\epsilon_{i}}(|y|_{+}^{2})>\epsilon\implies
M_{p_{i},w_{i}}^{\epsilon_{i}}(|y|_{-}^{2})=\epsilon$,

  * •

$M_{p_{i},w_{i}}^{\epsilon_{i}}(|y|_{-}^{2})>\epsilon\implies
M_{0,w_{i}}^{\epsilon_{i}}(|y|_{+}^{2})=\epsilon$.

Therefore,

$\displaystyle\min(y)<0$ $\displaystyle\iff\exists j\in[1{:}n]\text{ such that
}y_{j}<0$ $\displaystyle\iff
M_{p_{i},w_{i}}^{\epsilon_{i}}(|y|_{-}^{2})>\epsilon$ $\displaystyle\iff
M_{0,w_{i}}^{\epsilon_{i}}(|y|_{+}^{2})-M_{p_{i},w_{i}}^{\epsilon_{i}}(|y|_{-}^{2})<0$
$\displaystyle\iff{}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}}(y)<0$
$\displaystyle\min(y)>0$ $\displaystyle\iff y_{j}>0\text{ for all
}j\in[1{:}n]$ $\displaystyle\iff
M_{0,w_{i}}^{\epsilon_{i}}(|y|_{+}^{2})>\epsilon$ $\displaystyle\iff
M_{0,w_{i}}^{\epsilon_{i}}(|y|_{+}^{2})-M_{p_{i},w_{i}}^{\epsilon_{i}}(|y|_{-}^{2})>0$
$\displaystyle\iff{}^{\wedge}h_{p_{i},w_{i}}^{\epsilon_{i}}(y)>0$

######  Lemma 3.6.

(Limit behavior)

  * •

$\lim_{p\to\infty}\;{}^{\wedge}h_{p,w}^{\epsilon}(y)=\min(y)+\epsilon^{0.5}$
if $\min(y)<0$ and $\epsilon^{0.5}<|\min(y)|$.

  * •

$\lim_{p\to\infty}\;{}^{\vee}h_{p,w}^{\epsilon}(y)=\max(y)-\epsilon^{0.5}$ if
$\max(y)>0$ and $\epsilon^{0.5}<\max(y)$.

Proof: If $z\in\mathbb{R}_{+}$, then $M_{\infty,w}^{0}(z)=\max(z)$ for any
$w\in\mathbb{R}_{++}^{n}$ [43]. If $\epsilon^{0.5}<\max(z)$, then
$M_{\infty,w}^{\epsilon}(z)=\max(z)$. Therefore, if $\min(y)<0$ and
$\epsilon^{0.5}<|\min(y)|$, then
$\lim_{p\to\infty}\;{}^{\wedge}h_{p,w}^{\epsilon}(y)=\epsilon^{0.5}-\lim_{p\to\infty}\Big{(}M_{p,w}^{\epsilon}(|y|_{-}^{2})\Big{)}^{1/2}=\epsilon^{0.5}-(\max(|y|_{-}^{2}))^{1/2}=\epsilon^{0.5}-|\max(|y_{i}|_{-})|=\min(y)+\epsilon^{0.5}$.

Since ${}^{\vee}h_{p,w}^{\epsilon}(y):=-{}^{\wedge}h_{p,w}^{\epsilon}(-y)$,
$\lim_{p\to\infty}\;{}^{\vee}h_{p,w}^{\epsilon}(y)=-\min(-y)-\epsilon^{0.5}=\max(y)-\epsilon^{0.5}$
if $\max(y)>0$ and $\epsilon^{0.5}<\max(y)$.

###  3.2 Gradient properties of the D-GMSR

As discussed in Section 2.3, the robustness value of any STL specification is
determined by the value of a specific predicate function at a single time
instant in D-SR since all predicate functions included in the STL
specification and their values at all times are composed with the $\min$ and
$\max$ functions. Accordingly, when maximizing the robustness value of a
specification using D-SSR and a gradient-based optimization method, typically
only the gradient of a specific predicate function at a single time instant is
considered notably in each iteration. Conversely, key functions of D-GMSR,
denoted as ${}^{\wedge}h_{p,w}^{\epsilon}$ and ${}^{\vee}h_{p,w}^{\epsilon}$,
are formed through a specific composition of generalized means in (3). This
construction enables the simultaneous consideration of values from multiple
predicate functions across multiple time instances, addressing the issues of
locality and masking observed in D-SR. Accordingly, the gradient of an STL
specification is the weighted summation of the gradient of multiple predicate
functions at multiple time instances. These characteristics of key functions
of D-GMSR lead to the alleviation of certain undesirable optimization results
observed in D-SSR, where the STL specification is not met, as demonstrated by
an example in Section 4.

To demonstrate the characteristic of key functions of D-GMSR, suppose that
$f_{i}:\mathbb{R}^{n}\to\mathbb{R}$ is a function of $x\in\mathbb{R}^{n}$ for
all $i\in[1{:}m]$ and $y=(f_{1}(x),f_{2}(x),\dots,f_{m}(x))$, then

$\displaystyle\nabla_{x}{}^{\wedge}h_{p,w}^{\epsilon}($ $\displaystyle
y)=c_{0}^{+}\nabla_{x}M_{0,w}^{\epsilon}(|y|_{+}^{2})-c_{p}^{-}\nabla_{x}M_{p,w}^{\epsilon}(|y|_{-}^{2}),$
$\displaystyle\nabla_{x}{}^{\vee}h_{p,w}^{\epsilon}($ $\displaystyle
y)=c_{p}^{+}\nabla_{x}M_{p,w}^{\epsilon}(|y|_{+}^{2})-c_{0}^{-}\nabla_{x}M_{0,w}^{\epsilon}(|y|_{-}^{2}),$

where

$\displaystyle\nabla_{x}M_{p,w}^{\epsilon}(|y|_{\square}^{2})$
$\displaystyle=c_{p,m}^{\square}\sum_{i=1}^{m}w_{i}|f_{i}(x)|_{\square}^{2p-1}\nabla_{x}f_{i}(x),$
$\displaystyle\nabla_{x}M_{0,w}^{\epsilon}(|y|_{\square}^{2})$
$\displaystyle=c_{0,m}^{\square}\sum_{i=1}^{m}w_{i}\frac{1}{|f_{i}(x)|_{\square}}\nabla_{x}f_{i}(x),$
$\displaystyle c_{0}^{\square}$
$\displaystyle=\frac{1}{2}\Big{(}M_{0,w}^{\epsilon}(|y|_{\square}^{2})\Big{)}^{-1/2},$
$\displaystyle c_{p}^{\square}$
$\displaystyle=\frac{1}{2}\Big{(}M_{p,w}^{\epsilon}(|y|_{{\square}}^{2})\Big{)}^{-1/2},$
$\displaystyle c_{0,m}^{\square}$
$\displaystyle=\frac{2\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\prod_{j=1}^{m}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}|f_{j}(x)|_{\square}^{2w_{i}}}{\bm{1}^{T}w\Big{(}M_{0,w}^{\epsilon}(|y|_{\square}^{2})\Big{)}^{\bm{1}^{T}w-1}},$
$\displaystyle c_{p,m}^{\square}$
$\displaystyle=\frac{2p}{p\bm{1}^{T}w\Big{(}M_{p,w}^{\epsilon}(|y|_{\square}^{2})\Big{)}^{p-1}}.$

Consequently,

$\displaystyle\nabla_{x}{}^{\wedge}h_{p,w}^{\epsilon}(y)$
$\displaystyle=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla_{x}{}^{\wedge}h_{p,w}^{\epsilon}((f_{1}(x),f_{2}(x),\dots,f_{m}(x)))\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
$\displaystyle=\sum_{i=1}^{m}w_{i}a_{i}\nabla_{x}f_{i}(x)\text{ and}$
$\displaystyle\nabla_{x}{}^{\vee}h_{p,w}^{\epsilon}(y)$
$\displaystyle=\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\nabla_{x}{}^{\vee}h_{p,w}^{\epsilon}((f_{1}(x),f_{2}(x),\dots,f_{m}(x)))\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
$\displaystyle=\sum_{i=1}^{m}w_{i}b_{i}\nabla_{x}f_{i}(x)$

where

$\displaystyle a_{i}$
$\displaystyle=\begin{cases}c_{0}^{+}c_{0,m}^{+}\displaystyle\frac{1}{|f_{i}(x)|_{+}}&\text{if
}{}^{\wedge}h_{p,w}^{\epsilon}(y)>0\\\\[8.5359pt]
c_{p}^{-}c_{p,m}^{-}-|f_{i}(x)|^{2p-1}_{-}&\text{otherwise}\end{cases},$
$\displaystyle b_{i}$
$\displaystyle=\begin{cases}c_{0}^{-}c_{0,m}^{-}\displaystyle\frac{1}{-|f_{i}(x)|_{-}}&\text{if
}{}^{\vee}h_{p,w}^{\epsilon}(y)<0\\\\[8.5359pt]
c_{p}^{+}c_{p,m}^{+}|f_{i}(x)|^{2p-1}_{+}&\text{otherwise}\end{cases}.$

Therefore,

  * •

If ${}^{\wedge}h_{p,w}^{\epsilon}(y)>0$, then
$\frac{a_{i}}{a_{j}}\propto\frac{f_{j}(x)}{f_{i}(x)}$ for all $i,j\in[1{:}m]$.

  * •

If ${}^{\wedge}h_{p,w}^{\epsilon}(y)<0$, then
$\frac{a_{i}}{a_{j}}\propto\frac{|f_{i}(x)|^{2p-1}}{|f_{j}(x)|^{2p-1}}$ for
all $i,j\in\\{k:f_{k}<0,k\in[1{:}m]\\}$, and $a_{i}=0$ for all
$i\in\\{k:f_{k}\geq 0,k\in[1{:}m]\\}$.

  * •

If ${}^{\vee}h_{p,w}^{\epsilon}(y)<0$, then
$\frac{b_{i}}{b_{j}}\propto\frac{|f_{j}(x)|}{|f_{i}(x)|}$ for all
$i,j\in[1{:}m]$.

  * •

If ${}^{\vee}h_{p,w}^{\epsilon}(y)>0$, then
$\frac{b_{i}}{b_{j}}\propto\frac{f_{i}(x)^{2p-1}}{f_{j}(x)^{2p-1}}$ for all
$i,j\in\\{k:f_{k}>0,k\in[1{:}m]\\}$, and $b_{i}=0$ for all $i\in\\{k:f_{k}\leq
0,k\in[1{:}m]\\}$.

Recall that for the $\min$ and ${}^{\wedge}h_{p,w}^{\epsilon}$ functions to be
positive, all their variables must be positive. The derivative of the
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
function yields almost zero for all the variables that do not have the minimum
value. On the other hand, if the values of all the variables are negative,
then the derivative of the ${}^{\wedge}h_{p,w}^{\epsilon}$ function with
respect to each variable is non-zero and it increases as the value of the
variable decreases. In other words, the variable that is farthest from being
positive receives greater consideration during the optimization process while
other negative-valued variables are still not ignored. If there exist both
positive and negative-valued variables, then the derivative with respect to
negative-valued variables increases as the value of the variable decreases and
the derivative of the function with respect to non-negative-valued variables
becomes zero. In other words, only the variables causing the negative yield of
the ${}^{\wedge}h_{p,w}^{\epsilon}$ function are taken into account.

For the $\max$ and ${}^{\vee}h_{p,w}^{\epsilon}$ functions to be positive,
only one of their variables must be positive. The derivative of the
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
function yields almost zero for all the variables that do not have the maximum
value. On the other hand, if the values of all the variables are negative,
then the derivative of the ${}^{\vee}h_{p,w}^{\epsilon}$ function with respect
to each variable is non-zero and it increases as the value of the variable
increases. In other words, the variable that is closest to being positive
receives greater consideration during the optimization process while other
negative-valued variables are still not ignored. If there exist both positive
and negative-valued variables, then the derivative with respect to positive-
valued variables increases as the value of the variable increases and the
derivative of the function with respect to non-positive-valued variables
becomes zero. In other words, only the variables causing the positive yield of
the ${}^{\wedge}h_{p,w}^{\epsilon}$ function are taken into account.

As a result, unlike the
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
and
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
functions, while the values of the ${}^{\wedge}h_{p,w}^{\epsilon}$ and
${}^{\vee}h_{p,w}^{\epsilon}$ functions are negative, their derivatives with
respect to the variables whose variations can make these function positive are
never zero. Hence, certain undesirable optimization results observed in D-SSR,
where the STL specification is not met, are mitigated, as demonstrated by an
example in Section 4.

######  Example 3.7.

Similar to the example presented in Figure 1, we illustrate the variations in
the predicate function values during the maximization of
$\Gamma^{\varphi_{1}\wedge\varphi_{2}\wedge\dots\wedge\varphi_{5}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)$
and
$\Gamma^{\varphi_{1}\vee\varphi_{2}\vee\dots\vee\varphi_{5}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)$
functions in Figure 2 to provide numerical evidence for the statements
regarding the derivative of ${}^{\wedge}h_{p,w}^{\epsilon}$ and
${}^{\vee}h_{p,w}^{\epsilon}$ functions for $\epsilon=10^{-8}$, $p=1$ and
$w=\bm{1}$. Here, $\Gamma^{\varphi_{i}}(x,k)=f_{i}(x_{k})\in\mathbb{R}$ and
$f_{i}:\mathbb{R}\to\mathbb{R}$ for all $i\in[1{:}5]$. Therefore,
$\Gamma^{\varphi_{1}\wedge\varphi_{2}\wedge\dots\wedge\varphi_{5}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)={}^{\wedge}h_{p,w}^{\epsilon}(y)$
and
$\Gamma^{\varphi_{1}\vee\varphi_{2}\vee\dots\vee\varphi_{5}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)={}^{\vee}h_{p,w}^{\epsilon}(y)$,
where $y=(f_{1}(x_{k}),f_{2}(x_{k}),\dots,f_{5}(x_{k}))$. Alternatively,
maximizing $\Gamma^{\bm{G}_{[1:5]}\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,0)$
and $\Gamma^{\bm{F}_{[1:5]}\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,0)$
functions, respectively, where $\Gamma^{\varphi}(x,k)=f(x_{k})\in\mathbb{R}$
for all $k\in[1{:}5]$ and $f:\mathbb{R}\to\mathbb{R}$, leads to similar
changes in predicate function values across different time instances because
$\Gamma^{\bm{G}_{[1:5]}\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,0)={}^{\wedge}h_{p,w}^{\epsilon}(y)$
and
$\Gamma^{\bm{F}_{[1:5]}\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,0)={}^{\vee}h_{p,w}^{\epsilon}(y)$,
where $y=(f(x_{1}),f(x_{2}),\dots,f(x_{5}))$.

For ${}^{\wedge}h_{p,w}^{\epsilon}$ function, when there exist both positive
and negative-valued variables, the derivative of the function with respect to
its variables causes only the value of the negative-valued variables to
increase such that their value will become equal to each other at $0$,
regardless of their initial values. When there exists only positive-valued
variables, all variables increase such that their values will become equal to
each other since the variable that has a smaller value increases more. For
${}^{\vee}h_{p,w}^{\epsilon}$ function, when there exist only negative-valued
variables, all of them increase such that the variable that is closest to
being positive increases more. When there exist both positive and negative-
valued variables, only positive-valued variables increase such that the
variable with higher value increases more than other positive-valued
variables. These results validate the statements regarding the derivatives of
${}^{\wedge}h_{p,w}^{\epsilon}$ and ${}^{\vee}h_{p,w}^{\epsilon}$ functions.

Figure 2:  The changes in the values of $f_{i}(x_{k})$ (or $f(x_{i})$) and the
corresponding robustness values while maximizing the conjunction (or always)
and disjunction (or eventually) operators using D-GMSR are presented on the
left and right sides of the figure, respectively. The color gradient, ranging
from purple to red in the figures above, represents the values of the
functions from $f_{1}(x_{k})$ to $f_{5}(x_{k})$ (or $f(x_{1})$ to $f(x_{5})$),
respectively.

##  4 Numerical Simulations

Sequential convex programming (SCP) is a prominent direct method to solve
trajectory optimization problems [38, 39, 40, 41]. SCP algorithms are also
utilized to initialize indirect methods for trajectory optimization [44, 45].
The real-time capability of this approach is demonstrated for powered descent
guidance problems in [46, 47]. SCP algorithms also demonstrate good
performance in solving optimal control problems with state-triggered
constraints, which are a subset of STL constraints [48, 49, 50].

For solving the trajectory optimization problems we present, we use the
successive convexification framework [41], which combines multiple-shooting
for the exact discretization of the continuous-time dynamics [51, 52],
$\ell_{1}$ exact penalization for the nonconvex constraints [53, Chap. 17] and
a convergence-guaranteed SCP algorithm, prox-linear method [54], for solving
the resulting optimization problem. For modeling and solving each convex
subproblem encountered in every iteration of the prox-linear algorithm, we use
CVXPY [55] with either ECOS [56] or MOSEK [57]. Customizable conic solvers
such as the proportional-integral projected gradient (PIPG) algorithm [58, 59]
could also be employed for efficient implementations. We employ the adaptive
weight update algorithm [60, Alg. 2.2] to dynamically determine the proximal
term weight within the prox-linear algorithm. The convergence guarantee of the
prox-linear method, achieved through an appropriate selection of the proximal
term weight, ensures that the converged solutions, which are feasible with
respect to the finite-dimensional nonconvex optimization problem, are also
Karush-Kuhn-Tucker (KKT) points [41, Thm. 19]. The numerical performance of
the successive convexification framework is demonstrated for applications
ranging from GPU-accelerated trajectory optimization for 6-DoF rocket landing
[61] and nonlinear model predictive control (NMPC) for obstacle avoidance
[62], to trajectory optimization for 6-DoF aircraft approach and landing [63].

######  Remark 3.

A reformulation technique is proposed for the successive convexification
framework in [41] to ensure the satisfaction of the path constraints in
continuous time; however, it covers only always specification in STL. The
application of this approach to STL specifications in general has not been
explored. Therefore, given that we solely considered STL specifications at
node points, we impose path constraints exclusively at node points for
consistency.

###  4.1 Locality & Masking

We first solve a simple trajectory optimization problem in (5) to illustrate
an undesirable optimization result, where the STL specification is not met due
to the properties of
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\min}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
and
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\max}_{\kappa}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$
functions in D-SSR, followed by a demonstration of D-GMSR’s performance on the
same problem. The position and control input of the vehicle at time $k$ are
represented as $x_{k}$ and $u_{k}$, respectively. The control input of the
vehicle is upper-bounded and there exist boundary conditions on the initial
and final positions of the vehicle. The objective is for the vehicle to
eventually reach a circular area centered at $p_{w}$ with a radius of $r_{w}$.
The corresponding predicate function and the specification are defined as
$f(x):=r_{w}-\|x-p_{w}\|$ and $\varphi:=(f(x)\geq 0)$, respectively. The
solutions of the optimization problems utilize D-SSR for $\kappa=25$ [31] and
D-GMSR for $\epsilon=10^{-8}$, $p=1$, $w=\bm{1}$ are presented in Fig. (3).

$\displaystyle\underset{x_{k},u_{k}}{\text{maximize}}$
$\displaystyle\tilde{\rho}^{\bm{F}_{[1{:}K]}\varphi}(x,0)\text{ or
}\Gamma^{\bm{F}_{[1{:}K]}\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,0)$ (5)
subject to $\displaystyle x_{k+1}=x_{k}+u_{k},\;\forall k\in[1{:}K-1]$
$\displaystyle x_{1}=(0,0)$ $\displaystyle x_{K}=(8,0)$
$\displaystyle\|u_{k}\|\leq 1.5,\;\forall k\in[1{:}K]$

The total number of nodes ($K$) is set as $9$ and the initial trajectory for
the optimization problem is set as the black dotted curve in Fig. 3. In order
to satisfy the STL specification, at least one of the node points ($x_{k}$)
must be inside of the blue circle in Fig. 3. For the initial trajectory, the
second node point is the closest node point to the center of the blue circle,
thus $\partial_{x_{k}}\tilde{\rho}^{\bm{F}_{[1{:}K]}\varphi}(x,0)\approx 0$
for all $k\in[1{:}K]\setminus\\{2\\}$. In other words, when the D-SSR is used,
the optimizer only forces the second node point to converge to the center of
the blue circle. However, due to the control-upper bound, the distance between
two node points cannot be greater than $1.5$. Therefore, $x_{2}$ cannot
converge to the inside of the blue circle since $x_{1}=(0,0)$. As a result,
the initial trajectory converges to the trajectory presented via the red
dotted curve in Fig. 3. Even though alternative node points could satisfy the
specification, D-SSR’s focus solely on the second node point prevents the
satisfaction of the specification.

In D-GMSR, however, for the initial trajectory, we have

$\partial_{x_{k}}\Gamma^{\bm{F}_{[1{:}K]}\varphi}_{\bm{\epsilon},\bm{p},\bm{w}}(x,0)=a_{k}\partial_{x_{k}}f(x)$

where $a_{2}>a_{k}\gg 0$ for all $k\in[1{:}K]\setminus\\{2\\}$. In other
words, when D-GMSR is used, the optimizer prioritizes the convergence of the
second node point towards the center of the blue circle but it also enforces
other node points inversely proportional to their distances from the circular
area. Therefore, despite the control upper bound preventing the convergence of
the second node point within the blue circle, the optimization process ensures
that other node points continue to converge towards the circle’s center,
eventually resulting in one of them reaching the center. As a result, when
D-GMSR is used, the initial trajectory converges to the trajectory presented
via the blue dotted curve in Fig. 3.

Figure 3:  The initial and converged trajectories obtained using D-SSR and
D-GMSR, along with the changes in predicate function for each node point, and
the corresponding robustness values throughout the iterations, are
illustrated. In the position plot, the blue circle shows the area that needs
to be visited. The black dotted curve represents the initial trajectory, while
the red and blue dotted curves show the converged trajectories obtained using
D-SSR and D-GMSR, respectively. Node points of each trajectory, from $x_{1}$
to $x_{9}$, are represented with varying colors ranging from purple to red.
Node points are connected by dotted curves for clarity. Similarly, in the
central figures, the color gradient, ranging from purple to red, represents
the value of the predicate function for each node point from $f(x_{1})$ to
$f(x_{9})$, respectively.

###  4.2 Quadrotor flight

We address an optimal control problem for a point mass model of a quadrotor,
incorporating STL specifications. In this scenario, the objective for the
quadrotor is to fly to a designated destination while satisfying the
constraints on speed, tilt angle, and thrust force. However, the vehicle’s
speed must not exceed a specified threshold until it stops at the battery
charging station for a consecutive number of time steps.

####  4.2.1 Vehicle dynamics

The state and control input of the vehicle are defined as $x:=(r,v)$ and
$u:=T$, where $r:=(r^{x},r^{y},r^{z})\in\mathbb{R}^{3}$,
$v:=(v^{x},v^{y},v^{z})\in\mathbb{R}^{3}$, and
$T:=(T^{x},T^{y},T^{z})\in\mathbb{R}^{3}$ are the position, velocity, and
thrust force, respectively. The dynamic equation of the vehicle is defined as

$\displaystyle\dot{x}(t)$ $\displaystyle=f(x(t),u(t))$ (6)
$\displaystyle=\begin{bmatrix}v(t)\\\
T(t)/m-c_{\mathrm{d}}\|v(t)\|v(t)+g\end{bmatrix}$

where $m\in\mathbb{R}_{++}$ and $c_{\mathrm{d}}\in\mathbb{R}_{++}$ are the
mass and the drag coefficient of the vehicle, respectively, and
$g\in\mathbb{R}^{3}$ is the gravity vector.

We discretize the continuous-time dynamics using the multiple-shooting
approach. The time is discretized over a finite grid as
$0=t_{1}<t_{2}<\dots<t_{K}=t_{f}$, $x_{k}$ and $u_{k}$ represent the state and
control input at the node point $t_{k}$ for $k\in[1{:}K]$. The control input
$u$ is parameterized via first-order hold and the continuous-time dynamics (6)
is discretized using multiple-shooting as follows:

$\displaystyle x_{k+1}=F_{k}(x_{k},u_{k},u_{k+1})$ (7)

for $k\in[1{:}K-1]$. See [41, 62] for further details.

The boundary conditions are specified as follows:

$\displaystyle
r_{1}=r_{\mathrm{i}},\;r_{K}=r_{\mathrm{f}},\;v_{1}=v_{\mathrm{i}},\;v_{K}=v_{\mathrm{f}}$
(8a-d)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!u_{1}=mg,\;u_{K}=mg$
(8e-f)

The speed, the tilt angle, and the thrust force of the vehicle are constrained
as follows:

$\displaystyle\|v_{k}\|_{2}$ $\displaystyle\leq v_{\mathrm{max}}\;\;\;\forall
k\in[1{:}K]$ (9a) $\displaystyle\|u_{k}\|_{2}\cos(\theta)$ $\displaystyle\leq
e_{1}^{\top}u_{k}\;\;\;\forall k\in[1{:}K]$ (9b) $\displaystyle\|T_{k}\|_{2}$
$\displaystyle\leq T_{\mathrm{max}}\;\;\;\forall k\in[1{:}K]$ (9c)

####  4.2.2 STL specifications

In this scenario, the vehicle’s speed must not exceed a specified threshold
$v_{\mathrm{save}}$ until it stops at the battery charging station, a circular
area centered at $r_{\mathrm{w}}$ with radius $d_{\mathrm{w}}$, for
$k_{\mathrm{w}}$ consecutive time steps.

In formulating the corresponding STL specification, we initially define two
predicate functions along with their respective predicates as follows:

$\displaystyle f(x):=v_{\mathrm{save}}-\|v\|,\;\;$
$\displaystyle\;\;\varphi_{1}:=(f(x)\geq 0)$ $\displaystyle
g(x):=d_{\mathrm{w}}-\|r-r_{w}\|,\;\;$
$\displaystyle\;\;\varphi_{2}:=(g(x)\geq 0)$

to assess whether the vehicle’s speed exceeds the specified threshold and
whether the vehicle is positioned at the battery charging station,
respectively.

We subsequently define two predicates as follows

$\displaystyle\vartheta_{1}$
$\displaystyle:=\bm{G}_{[0{:}k_{\mathrm{w}}-1]}\varphi_{1}$
$\displaystyle\vartheta_{2}$
$\displaystyle:=\bm{G}_{[0{:}k_{\mathrm{w}}-1]}\varphi_{2}$

to assess if these specifications are satisfied for $k_{\mathrm{w}}$
consecutive time steps.

Given the scenario’s requirement to satisfy the speed specification, until the
vehicle is positioned at the battery charging station for $k_{\mathrm{w}}$
consecutive time steps, the resultant specification is formulated using D-GMSR
as follows:

$\Gamma^{\vartheta_{1}\bm{U}_{[0:K-(k_{\mathrm{w}}-1)]}\vartheta_{2}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,1)$
(10)

The construction of the specification using D-GMSR is presented in Appendix A.

The resulting problem is presented in (11) where Table 3 shows the system and
simulation parameters.

$\displaystyle\underset{x_{k},u_{k}}{\text{maximize}}$
$\displaystyle\eqref{eq:qf_stl}\quad\text{subject to}\quad$
$\displaystyle\eqref{eq:quad_dyn_dt},\eqref{eq:quad_bc},\eqref{eq:quad_cons}$
(11)

Table 3:  Parameter | Value
---|---
$t_{f}$, $K$ |  $20$ s, $21$
$r_{\mathrm{i}}$, $r_{\mathrm{f}}$ |  $(0,0,0)$ m, $(8,0,0)$ m
$v_{\mathrm{i}}$, $v_{\mathrm{f}}$ |  $0$ m s-1, $0$ m s-1
$T_{\mathrm{max}}$, $m$ |  $10.3$ N, $1.0$ kg
$v_{\mathrm{max}}$, $\theta_{\mathrm{max}}$ |  $2$ m s-1, $10$ deg
$v_{\mathrm{save}}$, $k_{\mathrm{w}}$ |  $1$ m s-1, 11
$r_{\mathrm{w}}$, $d_{\mathrm{w}}$ |  $(1.25,2,2)$ m, $0.2$ m
$g$, $c_{\mathrm{d}}$ |  $(0,0,-9.806)$ m s-2, $0.5$ m-1
$\epsilon_{i}\in\bm{\epsilon}$, $p_{i}\in\bm{p}$, $w_{i}\in\bm{w}$ |  $10^{-8}$, $1$, $\bm{1}$ for all $i\in\mathbb{Z}_{++}$

Since (7) is a nonconvex constraint, it is penalized with $\ell_{1}$ norm, and
the resulting optimization problem is solved via the prox-linear method. The
result, where the quadrotor executes the control input solution provided by
the optimization algorithm, is presented in Figure 4. While the constraints on
thrust, tilt angle, and speed are satisfied at the node points, inter-sample
constraint violations discussed in Remark 3 are observed on the speed
constraint. In conjunction with the reformulation technique proposed in [41],
continuous-time formulation of the STL specifications is required to ensure
the continuous-time satisfaction of both STL specifications and other
constraints. The quadrotor enters the battery charging station at the $5$th
time step and stays there until the $15$th time step, ensuring that its speed
remains below the threshold throughout this period. To maximize the robustness
value, the node points that are in the battery charging station converge
through the center of the station. Also, it is observed that the vehicle’s
speed remains significantly lower than the value of $v_{\textrm{save}}$ until
it exits the charging station.

Figure 4: Quadrotor flight with a battery charging station specification.

###  4.3 Autonomous Rocket landing

We address a powered descent guidance problem, specifically, the rocket
landing scenario presented in [64], by incorporating implication
specifications, i.e., state-triggered constraints (STCs). In this scenario,
the rocket must execute a precision landing, involving a maneuver to
transition from a horizontal to vertical orientation, and reduce the number of
operating engines from $3$ to $1$, similar to Starship [65]. While the
constraint for decreasing the number of operating engines is triggered by
speed, multiple constraints related to the landing’s safety are triggered by
altitude.

####  4.3.1 Vehicle dynamics

The state vector of the vehicle is defined as

$x=(m,r,v,\theta,\omega,\delta)\in\mathbb{R}^{8}$

where $m\in\mathbb{R}^{1}_{++}$, $r:=(r^{x},r^{z})\in\mathbb{R}^{2}$,
$v:=(v^{x},v^{z})\in\mathbb{R}^{2}$, $\theta\in\mathbb{R}^{1}$,
$\omega\in\mathbb{R}^{1}$ and $\delta\in\mathbb{R}^{1}$ represent the mass,
position, velocity, tilt angle, angular velocity of the vehicle, and gimbal
angle of the vehicle’s engine, respectively. The control input of the vehicle
is defined as

$\displaystyle u=(T,\overset{\text{{\large{.}}}}{\delta})\in\mathbb{R}^{4}$

where $T:=(T^{1},T^{2},T^{3})\in\mathbb{R}_{+}^{3}$, with
$T^{i}\in\mathbb{R}_{+}$ representing the thrust magnitude of the $i$th
engine, and $\overset{\text{{\large{.}}}}{\delta}\in\mathbb{R}^{1}$denotes the
gimbal rate of the vehicle’s engine. The dynamic equation of the vehicle is
defined as

$\displaystyle\dot{x}(t)$ $\displaystyle=f(x(t),u(t))$ (12)
$\displaystyle=\begin{bmatrix}-\alpha\bm{1}^{T}T(t)\\\ v\\\
g_{\mathcal{I}}+\frac{1}{m(t)}(F_{\mathcal{I}}(t)+D_{\mathcal{I}}(t))\\\
\omega\\\
\frac{1}{J_{\mathcal{B}}(t)}(-l_{\mathrm{cm}}\bm{1}^{T}T(t)\sin(\delta(t))-l_{\mathrm{cp}}e_{2}^{\top}D_{\mathcal{B}}(t))\\\
\overset{\text{{\large{.}}}}{\delta}\end{bmatrix}$

where

$\displaystyle\alpha$ $\displaystyle:=\frac{1}{I_{\mathrm{sp}}g_{0}}$
$\displaystyle g_{\mathcal{I}}$ $\displaystyle:=(0,g_{0})$ $\displaystyle
F_{\mathcal{I}}(t)$
$\displaystyle:=\bm{1}^{T}T(t)\begin{bmatrix}\cos(\theta(t)+\delta(t))\\\
-\sin(\theta(t)+\delta(t))\end{bmatrix}$ $\displaystyle D_{\mathcal{I}}(t)$
$\displaystyle:=\mathcal{R}_{\mathcal{I}\leftarrow\mathcal{B}}(t)D_{\mathcal{B}}(t)$
$\displaystyle\mathcal{R}_{\mathcal{I}\leftarrow\mathcal{B}}(t)$
$\displaystyle:=\begin{bmatrix}\cos(\theta(t))&\sin(\theta(t))\\\
-\sin(\theta(t))&\cos(\theta(t))\end{bmatrix}$ $\displaystyle
D_{\mathcal{B}}(t)$
$\displaystyle:=-\frac{1}{2}\rho_{\mathrm{air}}S_{\mathrm{area}}\|v(t)\|_{2}C_{\mathrm{aero}}\mathcal{R}_{\mathcal{I}\leftarrow\mathcal{B}}^{\top}(t)v(t)$
$\displaystyle J_{\mathcal{B}}(t)$
$\displaystyle:=m(t)\bigg{(}\frac{l_{\mathrm{r}}^{2}}{4}+\frac{l_{\mathrm{h}}^{2}}{12}\bigg{)}$

Here, $\alpha\in\mathbb{R}_{++}$ is the thrust-specific fuel consumption,
$I_{\mathrm{sp}}\in\mathbb{R}_{++}$ is the specific impulse of the rocket
engine, $g_{0}\in\mathbb{R}_{++}$ is the average gravitational acceleration of
the planet, $F_{\mathcal{I}}(t)$ and $D_{\mathcal{I}}(t)$ are the thrust and
drag forces on the vehicle with respect to the inertial frame, respectively,
$\mathcal{R}_{\mathcal{I}\leftarrow\mathcal{B}}(t)\in\mathrm{SO}(2)$ is the
rotation matrix that maps coordinates in the body frame to the inertial frame,
$D_{\mathcal{B}}(t)$ is the drag force on the vehicle with respect to the body
frame, $\rho_{\mathrm{air}}\in\mathbb{R}_{++}$ is the atmospheric density,
$S_{\mathrm{area}}\in\mathbb{R}_{++}$ is the reference area,
$C_{\mathrm{aero}}:=\mathrm{diag}\\{c_{x},c_{z}\\}$ is the aerodynamic
coefficient matrix, where $c_{x}$, $c_{z}\in\mathbb{R}_{++}$ are the
aerodynamic coefficients along the body $x$ and $z$ axes, respectively,
$J_{\mathcal{B}}(t)\in\mathbb{R}_{++}$ is the moment of inertia of the vehicle
about the body $y$ axis, $l_{\mathrm{r}}\in\mathbb{R}_{++}$ and
$l_{\mathrm{h}}\in\mathbb{R}_{++}$ are the radius and height of the fuselage,
respectively, $l_{\mathrm{cm}}\in\mathbb{R}_{++}$ is the thrust moment-arm
(distance between the vehicle mass-center and the engine gimbal hinge point),
and $l_{\mathrm{cp}}\in\mathbb{R}_{+}$ is the aerodynamic moment-arm (the
distance between the vehicle mass-center and the center of pressure).

We discretize the continuous-time dynamics using the multiple-shooting
approach. The time is discretized over a finite grid as
$0=t_{1}<t_{2}<\dots<t_{K}=t_{f}$, $x_{k}$ and $u_{k}$ represent the state and
control input at the node point $t_{k}$ for $k\in[1{:}K]$. The control input
$u$ is parameterized via first-order hold and the continuous-time dynamics (6)
is discretized using multiple-shooting as follows:

$\displaystyle x_{k+1}=F_{k}(x_{k},u_{k},u_{k+1})$ (13)

for $k\in[1{:}K-1]$. See [41, 62] for further details.

The boundary conditions are specified as follows:

$\displaystyle m_{1}=m_{\mathrm{i}},\;m_{K}\geq
m_{\mathrm{dry}},\;r_{1}=r_{\mathrm{i}},\;r_{K}=r_{\mathrm{f}}$ (14a-d)
$\displaystyle
v_{1}=v_{\mathrm{i}},\;v_{K}=v_{\mathrm{f}},\;\theta_{1}=\theta_{\mathrm{i}},\;\theta_{K}=\theta_{\mathrm{f}}$
(14e-h)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\omega_{1}=\omega_{\mathrm{i}},\;\omega_{K}=\omega_{\mathrm{f}},\;\delta_{1}=\delta_{\mathrm{i}}$
(14i-k)

where $m_{\mathrm{dry}}$ is the dry mass of the vehicle.

The altitude, gimbal angle, thrust, and gimbal rate are constrained as
follows:

$\displaystyle r_{\mathrm{min}}^{x}$ $\displaystyle\leq r^{x}_{k}$ (15a)
$\displaystyle\|\delta_{k}\|$ $\displaystyle\leq\delta_{\mathrm{max}}$ (15b)
$\displaystyle T_{\mathrm{min}}$ $\displaystyle\leq T_{k}^{1}\leq
T_{\mathrm{max}}$ (15c) $\displaystyle 0$ $\displaystyle\leq T_{k}^{2}\leq
T_{\mathrm{max}}$ (15d) $\displaystyle 0$ $\displaystyle\leq T_{k}^{3}\leq
T_{\mathrm{max}}$ (15e)
$\displaystyle\|\overset{\text{{\large{.}}}}{\delta_{k}}\|$
$\displaystyle\leq\overset{\text{{\large{.}}}}{\delta}_{\mathrm{max}}$ (15f)

for $k\in[1{:}K-1]$, where $r_{\mathrm{min}}^{x}$ is the minimum altitude.

####  4.3.2 STL specifications

We impose $14$ different implication specifications i.e. STC, where
constraints are triggered by the speed and altitude, as follows:

$\displaystyle\varphi_{i}:=$
$\displaystyle\;(\|v_{k^{\prime}}\|_{2}>v_{\mathrm{trig}})\\!\implies\\!(T_{\mathrm{min}}\leq
T_{k}^{i+1})\text{ for $i\\!=\\!1,2$}$ $\displaystyle\varphi_{i+2}:=$
$\displaystyle\;(\|v_{k^{\prime}}\|_{2}<v_{\mathrm{trig}})\\!\implies\\!(T_{k}^{i+1}=0)\text{
for $i\\!=\\!1,2$}$ $\displaystyle\varphi_{5}:=$
$\displaystyle\;(\|v_{k^{\prime}}\|_{2}<v_{\mathrm{trig}})\\!\implies\\!(\|v_{k}\|_{2}+\epsilon_{\mathrm{v}}\leq
v_{\mathrm{trig}})$ $\displaystyle\varphi_{6}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(\|v_{k}\|\leq
v_{\mathrm{stc}})$ $\displaystyle\varphi_{7}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(\theta_{k}\leq\theta_{\mathrm{stc}})$
$\displaystyle\varphi_{8}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(\theta_{k}\geq-\theta_{\mathrm{stc}})$
$\displaystyle\varphi_{9}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(\omega_{k}\leq\omega_{\mathrm{stc}})$
$\displaystyle\varphi_{10}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(\omega_{k}\geq-\omega_{\mathrm{stc}})$
$\displaystyle\varphi_{11}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(\delta_{k}\leq\delta_{\mathrm{stc}})$
$\displaystyle\varphi_{12}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(\delta_{k}\geq-\delta_{\mathrm{stc}})$
$\displaystyle\varphi_{13}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(r_{z_{k}}\leq
r^{x}_{k}\tan(\gamma_{\mathrm{stc}}))$ $\displaystyle\varphi_{14}:=$
$\displaystyle\;(r_{k^{\prime}}^{x}<r_{\mathrm{trig}}^{x})\\!\implies\\!(r_{z_{k}}\geq-r^{x}_{k}\tan(\gamma_{\mathrm{stc}}))$

$\forall k\in[k^{\prime}{:}K]$, where $v_{\mathrm{trig}}$ and
$r_{\mathrm{trig}}^{x}$ are the velocity and altitude that trigger the
constraints, respectively, $\gamma_{\mathrm{max}}$ is the maximum glideslope
angle, and $\epsilon_{\mathrm{v}}$ is an arbitrarily small positive number.

All specifications are formulated via D-GMSR and the following constraint is
imposed on the optimization problem:

$\Gamma_{\bm{\epsilon},\bm{p},\bm{w}}^{\varphi_{i}}(x,1)\geq 0\text{ for
}i=1,2,\dots,14$ (16)

The objective function of this problem is to minimize the fuel consumption of
the rocket. The resulting optimization problem is presented in (17) where
Table 4 shows the system and simulation parameters.

$\displaystyle\underset{x_{k},u_{k}}{\text{maximize}}\;\;m_{K}\;\;\text{subject
to}\;\eqref{eq:rl_dyn_dt},\eqref{eq:rl_bc},\eqref{eq:rl_cons},\eqref{eq:rl_stl}$
(17)

Figure 5: Rocket landing with state-triggered constraints. Table 4:  Parameter | Value
---|---
$t_{f}$, $K$ |  $20$ s, $21$
$I_{\mathrm{sp}}$, $g_{0}$ |  $330$ s, $9.806$ m s-1
$l_{\mathrm{r}}$, $l_{\mathrm{h}}$ |  $4.5$ m, $50$ m
$l_{\mathrm{cm}}$, $l_{\mathrm{cp}}$ |  $20$ m, $10$ m
$\rho_{\mathrm{air}}$, $S_{\mathrm{area}}$ |  $1.225$ kg m-3, $545$ m2
$c_{x}$, $c_{z}$ |  $0.0522$, $0.4068$
$m_{\mathrm{i}}$, $m_{\mathrm{dry}}$ |  $100000$ kg, $85000$ kg
$r_{\mathrm{i}}$, $r_{\mathrm{f}}$ |  $(500,100)$ m, $(0,0)$ m
$v_{\mathrm{i}}$, $v_{\mathrm{f}}$ |  $(-85,0)$ m s-1, $(0,0)$ m s-1
$\theta_{\mathrm{i}}$, $\theta_{\mathrm{f}}$ |  $90$ deg, $0$ deg
$\omega_{\mathrm{i}}$, $\omega_{\mathrm{f}}$, $\delta_{\mathrm{i}}$ |  $0$ deg s-1, $0$ deg s-1, $0$ deg
$r_{\mathrm{min}}^{x}$, $\delta_{\mathrm{max}}$ |  $0$ m $10$ deg
$T_{\mathrm{max}}$, $T_{\mathrm{min}}$, $\overset{\text{{\large{.}}}}{\delta}_{\mathrm{max}}$ |  $2200$ kN, $880$ kN, $15$ deg s-1
$v_{\mathrm{trig}}$, $r_{\mathrm{trig}}^{x}$, $v_{\mathrm{stc}}$ |  $35$ m s-1, $100$ m, $20$ m s-1
$\theta_{\mathrm{stc}}$, $\omega_{\mathrm{stc}}$, $\delta_{\mathrm{stc}}$, $\gamma_{\mathrm{stc}}$ |  $5$ deg, $2.5$ deg s-1, $1$ deg, $5$ deg
$\epsilon_{i}\in\bm{\epsilon}$, $p_{i}\in\bm{p}$, $w_{i}\in\bm{w}$ |  $10^{-8}$, $1$, $\bm{1}$ for all $i\in\mathbb{Z}_{++}$

Since (13) and (16) are nonconvex constraints, they are penalized with
$\ell_{1}$ norm, and the resulting optimization problem is solved via the
prox-linear method. The result, where the rocket executes the control input
solution provided by the optimization algorithm, is presented in Figure 5.
While the constraints are satisfied at the node points, inter-sample
constraint violations discussed in Remark 3 are observed on the gimbal angle
constraint. In conjunction with the reformulation technique proposed in [41],
continuous-time formulation of the STL specifications is required to ensure
the continuous-time satisfaction of both STL specifications and other
constrains. During precision landing, the vehicle performs a maneuver to
transition from horizontal to vertical orientation, akin to Starship [65], and
it satisfies all constraints triggered by speed and altitude.

##  5 Conclusion and Future work

This paper introduces a $\mathcal{C}^{1}$-smooth, sound, and complete
robustness measure for STL specifications. The robustness measure demonstrates
desirable gradient properties and therefore addresses the locality and masking
problems, which are critical for numerical optimization. The performance of
the robustness measure is demonstrated on real-world trajectory optimization
problems for quadrotor flight and autonomous rocket landing solved within an
SCP framework.

The continuous-time modeling of STL specifications and a successive
convexification based solution method to ensure the continuous-time
satisfaction of STL specifications will be the subject of future studies.
Additionally, potential future work may involve the design of optimization
algorithms that exploit the inherent structure of STL specifications to
achieve faster convergence.  The authors gratefully acknowledge Fabio Spada
and Abhinav G. Kamath for helpful discussions and their feedback on the
initial drafts of the paper.

## References

  * [1] Christel Baier and Joost-Pieter Katoen.  Principles of Model Checking.  MIT press, 2008.  URL: https://mitpress.mit.edu/9780262026499/principles-of-model-checking/.
  * [2] Yunus Emre Sahin, Rien Quirynen, and Stefano Di Cairano.  Autonomous vehicle decision-making and monitoring based on signal temporal logic and mixed-integer programming.  In 2020 American Control Conference (ACC), pages 454–459. IEEE, 2020.  URL: https://doi.org/10.23919/ACC45564.2020.9147917.
  * [3] Vince Kurtz and Hai Lin.  Trajectory optimization for high-dimensional nonlinear systems under STL specifications.  IEEE Control Systems Letters, 5(4):1429–1434, 2020.  URL: https://doi.org/10.1109/LCSYS.2020.3038640.
  * [4] Franck Djeumou, Zhe Xu, Murat Cubuktepe, and Ufuk Topcu.  Probabilistic control of heterogeneous swarms subject to graph temporal logic specifications: A decentralized and scalable approach.  IEEE Transactions on Automatic Control, 68(4):2245–2260, 2022.  URL: https://doi.org/10.1109/TAC.2022.3176797.
  * [5] Yunus Emre Sahin, Petter Nilsson, and Necmiye Ozay.  Multirobot coordination with counting temporal logics.  IEEE Transactions on Robotics, 36(4):1189–1206, 2019.  URL: https://doi.org/10.1109/TRO.2019.2957669.
  * [6] Ali Tevfik Buyukkocak, Derya Aksaray, and Yasin Yazıcıoğlu.  Planning of heterogeneous multi-agent systems under signal temporal logic specifications with integral predicates.  IEEE Robotics and Automation Letters, 6(2):1375–1382, 2021.  URL: https://doi.org/10.1109/LRA.2021.3057049.
  * [7] Yash Vardhan Pant, Houssam Abbas, Rhudii A Quaye, and Rahul Mangharam.  Fly-by-logic: Control of multi-drone fleets with temporal logic objectives.  In 2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS), pages 186–197. IEEE, 2018.  URL: https://doi.org/10.1109/ICCPS.2018.00026.
  * [8] Bariş Başpinar, Hamsa Balakrishnan, and Emre Koyuncu.  Mission planning and control of multi-aircraft systems with signal temporal logic specifications.  IEEE Access, 7:155941–155950, 2019.  URL: https://doi.org/10.1109/ACCESS.2019.2949426.
  * [9] Meiyi Ma, Ji Gao, Lu Feng, and John Stankovic.  STLnet: Signal temporal logic enforced multivariate recurrent neural networks.  Advances in Neural Information Processing Systems, 33:14604–14614, 2020.  URL: https://dl.acm.org/doi/abs/10.5555/3495724.3496948.
  * [10] Derya Aksaray, Austin Jones, Zhaodan Kong, Mac Schwager, and Calin Belta.  Q-learning for robust satisfaction of signal temporal logic specifications.  In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 6565–6570. IEEE, 2016.  URL: https://doi.org/10.1109/CDC.2016.7799279.
  * [11] Mohammed Alshiekh, Roderick Bloem, Rüdiger Ehlers, Bettina Könighofer, Scott Niekum, and Ufuk Topcu.  Safe reinforcement learning via shielding.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018.  URL: https://doi.org/10.1609/AAAI.V32I1.11797.
  * [12] Karen Leung, Nikos Aréchiga, and Marco Pavone.  Backpropagation through signal temporal logic specifications: Infusing logical structure into gradient-based methods.  The International Journal of Robotics Research, 42(6):356–370, 2023\.  URL: https://doi.org/10.1177/02783649221082115.
  * [13] Amir Pnueli.  The temporal logic of programs.  In 18th Annual Symposium on Foundations of Computer Science (SFCS 1977), pages 46–57. IEEE, 1977.  URL: https://doi.org/10.1109/SFCS.1977.32.
  * [14] Rajeev Alur and David L Dill.  A theory of timed automata.  Theoretical Computer Science, 126(2):183–235, 1994.  URL: https://doi.org/10.1016/0304-3975(94)90010-8.
  * [15] Tichakorn Wongpiromsarn, Ufuk Topcu, and Richard M Murray.  Receding horizon temporal logic planning.  IEEE Transactions on Automatic Control, 57(11):2817–2830, 2012\.  URL: https://doi.org/10.1109/TAC.2012.2195811.
  * [16] Tichakorn Wongpiromsarn, Ufuk Topcu, Necmiye Ozay, Huan Xu, and Richard M Murray.  TuLiP: a software toolbox for receding horizon temporal logic planning.  In Proceedings of the 14th International Conference on Hybrid Systems: Computation and Control, pages 313–314, 2011.  URL: https://doi.org/10.1145/1967701.1967747.
  * [17] Sertac Karaman, Ricardo G Sanfelice, and Emilio Frazzoli.  Optimal control of mixed logical dynamical systems with linear temporal logic specifications.  In 2008 47th IEEE Conference on Decision and Control, pages 2117–2122. IEEE, 2008.  URL: https://doi.org/10.1109/CDC.2008.4739370.
  * [18] Ron Koymans.  Specifying real-time properties with metric temporal logic.  Real-Time Systems, 2(4):255–299, 1990.  URL: https://doi.org/10.1007/BF01995674.
  * [19] Rajeev Alur, Tomás Feder, and Thomas A Henzinger.  The benefits of relaxing punctuality.  Journal of the ACM (JACM), 43(1):116–146, 1996.  URL: https://doi.org/10.1145/227595.227602.
  * [20] Oded Maler and Dejan Nickovic.  Monitoring temporal properties of continuous signals.  In International Symposium on Formal Techniques in Real-Time and Fault-Tolerant Systems, pages 152–166. Springer, 2004.  URL: https://doi.org/10.1007/978-3-540-30206-3_12.
  * [21] Georgios E Fainekos and George J Pappas.  Robustness of temporal logic specifications for continuous-time signals.  Theoretical Computer Science, 410(42):4262–4291, 2009.  URL: https://doi.org/10.1016/j.tcs.2009.06.021.
  * [22] Alexandre Donzé and Oded Maler.  Robust satisfaction of temporal logic over real-valued signals.  In International Conference on Formal Modeling and Analysis of Timed Systems, pages 92–106. Springer, 2010.  URL: https://doi.org/10.1007/978-3-642-15297-9_9.
  * [23] Lars Lindemann and Dimos V Dimarogonas.  Control barrier functions for signal temporal logic tasks.  IEEE Control Systems Letters, 3(1):96–101, 2018.  URL: https://doi.org/10.1109/LCSYS.2018.2853182.
  * [24] Calin Belta and Sadra Sadraddini.  Formal methods for control synthesis: An optimization perspective.  Annual Review of Control, Robotics, and Autonomous Systems, 2:115–140, 2019.  URL: https://doi.org/10.1146/annurev-control-053018-023717.
  * [25] Vasumathi Raman, Alexandre Donzé, Mehdi Maasoumy, Richard M Murray, Alberto Sangiovanni-Vincentelli, and Sanjit A Seshia.  Model predictive control with signal temporal logic specifications.  In 53rd IEEE Conference on Decision and Control, pages 81–87. IEEE, 2014.  URL: https://doi.org/10.1109/CDC.2014.7039363.
  * [26] Sadra Sadraddini and Calin Belta.  Robust temporal logic model predictive control.  In 2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 772–779. IEEE, 2015.  URL: https://doi.org/10.1109/ALLERTON.2015.7447084.
  * [27] Alëna Rodionova, Lars Lindemann, Manfred Morari, and George J Pappas.  Time-robust control for STL specifications.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 572–579. IEEE, 2021.  URL: https://doi.org/10.1109/CDC45484.2021.9683477.
  * [28] Alëna Rodionova, Lars Lindemann, Manfred Morari, and George J Pappas.  Combined left and right temporal robustness for control under STL specifications.  IEEE Control Systems Letters, 7:619–624, 2022.  URL: https://doi.org/10.1109/LCSYS.2022.3209928.
  * [29] Yash Vardhan Pant, Houssam Abbas, and Rahul Mangharam.  Smooth operator: Control using the smooth robustness of temporal logic.  In 2017 IEEE Conference on Control Technology and Applications (CCTA), pages 1235–1240. IEEE, 2017.  URL: https://doi.org/10.1109/CCTA.2017.8062628.
  * [30] Iman Haghighi, Noushin Mehdipour, Ezio Bartocci, and Calin Belta.  Control from signal temporal logic specifications with smooth cumulative quantitative semantics.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 4361–4366. IEEE, 2019.  URL: https://doi.org/10.1109/CDC40024.2019.9029429.
  * [31] Yann Gilpin, Vince Kurtz, and Hai Lin.  A smooth robustness measure of signal temporal logic for symbolic control.  IEEE Control Systems Letters, 5(1):241–246, 2020.  URL: https://doi.org/10.1109/LCSYS.2020.3001875.
  * [32] Yuanqi Mao, Behcet Acikmese, Pierre-Loic Garoche, and Alexandre Chapoutot.  Successive convexification for optimal control with signal temporal logic specifications.  In Proceedings of the 25th ACM International Conference on Hybrid Systems: Computation and Control, pages 1–7, 2022.  URL: https://doi.org/10.1145/3501710.3519518.
  * [33] Lars Lindemann and Dimos V Dimarogonas.  Robust control for signal temporal logic specifications using discrete average space robustness.  Automatica, 101:377–387, 2019.  URL: https://doi.org/10.1016/j.automatica.2018.12.022.
  * [34] Noushin Mehdipour, Cristian-Ioan Vasile, and Calin Belta.  Arithmetic-geometric mean robustness for control from signal temporal logic specifications.  In 2019 American Control Conference (ACC), pages 1690–1695. IEEE, 2019.  URL: https://doi.org/10.23919/ACC.2019.8814487.
  * [35] Noushin Mehdipour, Cristian-Ioan Vasile, and Calin Belta.  Specifying user preferences using weighted signal temporal logic.  IEEE Control Systems Letters, 5(6):2006–2011, 2020.  URL: https://doi.org/10.1109/LCSYS.2020.3047362.
  * [36] Noushin Mehdipour, Cristian-Ioan Vasile, and Calin Belta.  Average-based robustness for continuous-time signal temporal logic.  In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 5312–5317. IEEE, 2019.  URL: https://doi.org/10.1109/CDC40024.2019.9029989.
  * [37] Takumi Akazaki and Ichiro Hasuo.  Time robustness in MTL and expressivity in hybrid system falsification.  In International Conference on Computer Aided Verification, pages 356–374. Springer, 2015.  URL: https://doi.org/10.1007/978-3-319-21668-3_21.
  * [38] Yuanqi Mao, Michael Szmuk, and Behçet Açıkmeşe.  Successive convexification of non-convex optimal control problems and its convergence properties.  In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 3636–3641. IEEE, 2016.  URL: https://doi.org/10.1109/CDC.2016.7798816.
  * [39] Danylo Malyuta, Yue Yu, Purnanand Elango, and Behçet Açıkmeşe.  Advances in trajectory optimization for space vehicle control.  Annual Reviews in Control, 52:282–315, 2021.  URL: https://doi.org/10.1016/j.arcontrol.2021.04.013.
  * [40] Danylo Malyuta, Taylor P Reynolds, Michael Szmuk, Thomas Lew, Riccardo Bonalli, Marco Pavone, and Behçet Açıkmeşe.  Convex optimization for trajectory generation: A tutorial on generating dynamically feasible trajectories reliably and efficiently.  IEEE Control Systems Magazine, 42(5):40–113, 2022.  URL: https://doi.org/10.1109/MCS.2022.3187542.
  * [41] Purnanand Elango, Dayou Luo, Abhinav G. Kamath, Samet Uzun, Taewan Kim, and Behçet Açıkmeşe.  Successive convexification for trajectory optimization with continuous-time constraint satisfaction.  arXiv preprint arXiv:2404.16826, 2024.  URL: https://doi.org/10.48550/arXiv.2404.16826.
  * [42] Danylo Malyuta and Behçet Açıkmeşe.  Fast homotopy for spacecraft rendezvous trajectory optimization with discrete logic.  Journal of Guidance, Control, and Dynamics, 46(7):1262–1279, 2023\.  URL: https://doi.org/10.2514/1.G006295.
  * [43] Harold Shniad.  On the convexity of mean value functions.  Bulletin of the American Mathematical Society, 54:770–776, 1948\.  URL: https://doi.org/10.1090/S0002-9904-1948-09077-2.
  * [44] Gao Tang, Fanghua Jiang, and Junfeng Li.  Fuel-optimal low-thrust trajectory optimization using indirect method and successive convex programming.  IEEE Transactions on Aerospace and Electronic Systems, 54(4):2053–2066, 2018.  URL: https://doi.org/10.1109/TAES.2018.2803558.
  * [45] Fabio Spada, Marco Sagliano, and Francesco Topputo.  Direct–indirect hybrid strategy for optimal powered descent and landing.  Journal of Spacecraft and Rockets, 60(6):1787–1804, 2023.  URL: https://doi.org/10.2514/1.A35650.
  * [46] Purnanand Elango, Abhinav Kamath, Yue Yu, John M. Carson III, and Behcet Acikmese.  A customized first-order solver for real-time powered-descent guidance.  In AIAA SciTech 2022 Forum. American Institute of Aeronautics and Astronautics, January 2022.  URL: https://doi.org/10.2514/6.2022-0951.
  * [47] Abhinav G. Kamath, Purnanand Elango, Skye Mceowen, Yue Yu, John M. Carson III, Mehran Mesbahi, and Behcet Acikmese.  Customized real-time first-order methods for onboard dual quaternion-based 6-DoF powered-descent guidance.  In AIAA SciTech 2023 Forum. American Institute of Aeronautics and Astronautics, January 2023.  URL: https://doi.org/10.2514/6.2023-2003.
  * [48] Michael Szmuk, Danylo Malyuta, Taylor P Reynolds, Margaret Skye Mceowen, and Behçet Açikmeşe.  Real-time quad-rotor path planning using convex optimization and compound state-triggered constraints.  In 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 7666–7673. IEEE, 2019.  URL: https://doi.org/10.1109/IROS40897.2019.8967706.
  * [49] Michael Szmuk, Taylor P Reynolds, and Behçet Açıkmeşe.  Successive convexification for real-time six-degree-of-freedom powered descent guidance with state-triggered constraints.  Journal of Guidance, Control, and Dynamics, 43(8):1399–1413, 2020\.  URL: https://doi.org/10.2514/1.G004549.
  * [50] Taylor P Reynolds, Michael Szmuk, Danylo Malyuta, Mehran Mesbahi, Behçet Açıkmeşe, and John M Carson III.  Dual quaternion-based powered descent guidance with state-triggered constraints.  Journal of Guidance, Control, and Dynamics, 43(9):1584–1599, 2020\.  URL: https://doi.org/10.2514/1.G004536.
  * [51] H G Bock and K J Plitt.  A multiple shooting algorithm for direct solution of optimal control problems.  IFAC proceedings volumes, 17(2):1603–1608, July 1984.  URL: https://doi.org/10.1016/s1474-6670(17)61205-9.
  * [52] Rien Quirynen, Milan Vukov, and Moritz Diehl.  Multiple shooting in a microsecond.  In Multiple Shooting and Time Domain Decomposition Methods: MuS-TDD, Heidelberg, May 6-8, 2013, pages 183–201. Springer, 2015.  URL: https://doi.org/10.1007/978-3-319-23321-5_7.
  * [53] Jorge Nocedal and Stephen Wright.  Numerical Optimization.  Springer, 2 edition, July 2006.  URL: https://doi.org/10.1007/978-0-387-40065-5.
  * [54] Dmitriy Drusvyatskiy and Courtney Paquette.  Efficiency of minimizing compositions of convex functions and smooth maps.  Mathematical Programming, 178:503–558, 2019.  URL: https://doi.org/10.1007/s10107-018-1311-3.
  * [55] Steven Diamond and Stephen Boyd.  CVXPY: A python-embedded modeling language for convex optimization.  Journal of Machine Learning Research, 17(83):1–5, 2016.  URL: https://www.cvxpy.org/.
  * [56] Alexander Domahidi, Eric Chu, and Stephen Boyd.  ECOS: An SOCP solver for embedded systems.  In 2013 European control conference (ECC), pages 3071–3076. IEEE, 2013.  URL: https://doi.org/10.23919/ECC.2013.6669541.
  * [57] Mosek ApS.  Mosek optimization toolbox for matlab.  User’s Guide and Reference Manual, Version, 4(1), 2019.  URL: https://docs.mosek.com/10.1/toolbox.pdf.
  * [58] Yue Yu, Purnanand Elango, Ufuk Topcu, and Behçet Açıkmeşe.  Proportional–integral projected gradient method for conic optimization.  Automatica, 142:110359, 2022.  URL: https://doi.org/10.1016/j.automatica.2022.110359.
  * [59] Yue Yu, Purnanand Elango, Behçet Açıkmeşe, and Ufuk Topcu.  Extrapolated proportional-integral projected gradient method for conic optimization.  IEEE Control Systems Letters, 7:73–78, 2022.  URL: https://doi.org/10.1109/LCSYS.2022.3186647.
  * [60] Coralia Cartis, Nicholas IM Gould, and Philippe L Toint.  On the evaluation complexity of composite function minimization with applications to nonconvex nonlinear programming.  SIAM Journal on Optimization, 21(4):1721–1739, 2011.  URL: https://doi.org/10.1137/11082381X.
  * [61] Govind M Chari, Abhinav G Kamath, Purnanand Elango, and Behçet Açıkmeşe.  Fast monte carlo analysis for 6-DoF powered-descent guidance via GPU-accelerated sequential convex programming.  In AIAA SciTech 2024 Forum. AIAA, 2024.  URL: https://doi.org/10.2514/6.2024-1762.
  * [62] Samet Uzun, Purnanand Elango, Abhinav G. Kamath, Taewan Kim, and Behçet Açıkmeşe.  Successive convexification for nonlinear model predictive control with continuous-time constraint satisfaction.  In 8th IFAC Conference on Nonlinear Model Predictive Control, Kyoto, Japan, August 2024.  URL: https://doi.org/10.48550/arXiv.2405.00061.
  * [63] Taewan Kim, Abhinav G. Kamath, Niyousha Rahimi, Jasper Corleis, Behçet Açıkmeşe, and Mehran Mesbahi.  Approach and landing trajectory optimization for a 6-dof aircraft with a runway alignment constraint, 2024.
  * [64] Abhinav G. Kamath, Purnanand Elango, Yue Yu, Skye Mceowen, Govind M. Chari, John M. Carson, III, and Behçet Açıkmeşe.  Real-time sequential conic optimization for multi-phase rocket landing guidance.  IFAC-PapersOnLine, 56(2):3118–3125, 2023.  URL: https://doi.org/10.1016/j.ifacol.2023.10.1444.
  * [65] Austin DeSisto.  Starship and its belly flop maneuver, 2021.  URL: https://everydayastronaut.com/starships-belly-flop-maneuver.

## Appendix

##  Appendix A Construction of the specification in quadrotor flight problem

We first define ${}^{f}s^{1}_{k}$ and ${}^{g}s^{1}_{k}$ functions to check
whatever the speed is less than $v_{\mathrm{save}}$ and the quadrotor is in
the battery charging station, respectively, from time step $k$ to the time
step $k+(k_{\mathrm{w}}-1)$. Note that  $\begin{aligned}
f(x)=v_{\mathrm{save}}-\|v\|,\;&\;\varphi_{1}:=(f(x)\geq
0),\;\vartheta_{1}=\bm{G}_{[0{:}k_{\mathrm{w}}-1]}\varphi_{1}\\\
g(x)=d_{\mathrm{w}}-\|r-r_{w}\|,\;&\;\varphi_{2}:=(g(x)\geq
0),\;\vartheta_{2}=\bm{G}_{[0{:}k_{\mathrm{w}}-1]}\varphi_{2}\end{aligned}$

Then

$\displaystyle{}^{f}s^{1}_{k}(x):=$
$\displaystyle\Gamma^{\vartheta_{1}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)$
$\displaystyle=$
$\displaystyle{}^{\wedge}h_{p_{1},w_{1}}^{\epsilon_{1}}((f(x_{k}),f(x_{k+1}),\dots,f(x_{k+k_{\mathrm{w}}-1})))$
$\displaystyle{}^{g}s^{1}_{k}(x):=$
$\displaystyle\Gamma^{\vartheta_{2}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,k)$
$\displaystyle=$
$\displaystyle{}^{\wedge}h_{p_{2},w_{2}}^{\epsilon_{2}}((g(x_{k}),g(x_{k+1}),\dots,g(x_{k+k_{\mathrm{w}}-1})))$

$\forall k\in[1{:}K-(k_{\mathrm{w}}-1)]$.

We now define $s_{k}^{2}$ function such that the output is positive if and
only if speed is less than $v_{\mathrm{save}}$ until the time step
$k+(k_{\mathrm{w}}-1)$ as follows:

$\displaystyle s_{k}^{2}(x)$
$\displaystyle:={}^{\wedge}h_{p_{3},w_{3}}^{\epsilon_{3}}(({}^{f}s^{1}_{1}(x),{}^{f}s^{1}_{2}(x),\dots,{}^{f}s^{1}_{k}(x)))$

$\forall k\in[1{:}K-(k_{\mathrm{w}}-1)]$.

We then define $s^{3}_{k}$ function such that the output is positive if and
only if speed is less than $v_{\mathrm{save}}$ until the time step
$k+(k_{\mathrm{w}}-1)$ and the being in a battery charging station
specification is satisfied between the time steps $k$ and
$k+(k_{\mathrm{w}}-1)$ as follows:

$\displaystyle s_{k}^{3}(x)$
$\displaystyle:={}^{\wedge}h_{p_{4},w_{4}}^{\epsilon_{4}}((s_{k}^{2}(x),{}^{g}s^{1}_{k}(x)))$

$\forall k\in[1{:}K-(k_{\mathrm{w}}-1)]$.

The $\vartheta_{1}\bm{U}_{[0:K-(k_{\mathrm{w}}-1)]}\vartheta_{2}$
specification is satisfied if and only if $s_{k}^{3}(x)>0$ for any
$k\in[1{:}K-(k_{\mathrm{w}}-1)]$; therefore,

$\displaystyle\Gamma^{\vartheta_{1}\bm{U}_{[0:K-(k_{\mathrm{w}}-1)]}\vartheta_{2}}_{\bm{\epsilon},\bm{p},\bm{w}}(x,1)={}^{\vee}h_{p_{5},w_{5}}^{\epsilon_{5}}(($
$\displaystyle s^{3}_{1}(x),s^{3}_{2}(x),\dots,$ $\displaystyle
s^{3}_{K-(k_{\mathrm{w}}-1)}(x)))$