Efficient Norm-Based Reachable Sets via
Iterative Dynamic Programming

Abstract

In this work, we present a numerical optimal control framework for reachable set computation using normotopes, a new set representation as a norm ball with a shaping matrix. In reachable set computations, we expect to continuously vary the shape matrix as a function of time. Incorporating the shape dynamics as an input, we build a controlled embedding system using a linear differential inclusion overapproximating the dynamics of the system, where a single forward simulation of this embedding system always provides an overapproximating reachable set of the system, no matter the choice of hypercontrol. By iteratively solving a linear quadratic approximation of the nonlinear optimal hypercontrol problem, we synthesize less conservative final reachable sets, providing a natural tradeoff between runtime and accuracy. Terminating our algorithm at any point always returns a valid reachable set overapproximation.

I Introduction

One way to assess the safety of a complex nonlinear system under uncertainty is to verify the specification for its forward reachable set. For nonlinear systems, the reachable set is often too complex to determine in closed form—instead, various algorithmic methods have been developed to compute guaranteed overapproximations. To combat conservativeness in these algorithmic strategies, different set representations have been developed, which have different tradeoffs between reachable set accuracy and computational efficiency. To name a few, there are zonotopes [1], constrained zonotopes [2], polynomial zonotopes [3], hybrid zonotopes [4], ellipsotopes [5], and taylor models [6]. For linear systems, any generalized star set [7] can be propagated using only $n+1$ simulations. Beyond these generator-based approaches, there are many other frameworks for reachable set computation. For instance, level set methods use the Hamilton-Jacobi equations [8, 9, 10] to represent the reachable set as a level set of a solution to a partial differential equation.

From a dual perspective, one can evolve halfspaces in linear systems using the adjoint equation [11]. This dual formulation was recently re-explored in [12], where they consider general parametric set representations defined by the intersection of constraints, instead of ones generated by constrained combinations. As shown in [12], the parameters defining the constraints can be arbitrarily updated as long as the offset term compensates accordingly. This observation is the basic backbone of this work—we directly optimize over the dynamics of these parameters to attempt to minimize the final volume of our overapproximating reachable set.

Computing reachable sets using norm balls is not a novel approach on its own, and there are many works applying results from contraction theory to compute such overapproximations. For instance, [13, 14] use bounds on logarithmic norms to bloat or shrink norm balls along nominal trajectories. The works [15, 16] use interval analysis to overapproximate the Jacobian matrix to automatically bound the logarithmic norm. More recently, [17] uses a more accurate interval matrix to bound the logarithmic norm by comparing specifically to the nominal trajectory. The advantage of these works is in the simplicity of computing the logarithmic norm, which is often known in closed form—for example, as a simple eigenvalue maximization in the $\ell_{2}$ case, or a maximum row/column sum computation in the $\ell_{\infty}$/$\ell_{1}$ cases [18, 19]. However, to our knowledge, existing works fix the norm along a real portion of the trajectory, using semi-definite programming to find a single shaping matrix along that segment. In this work, instead of successively choosing fixed norms for segments of the trajectory using semi-definite programming, we consider the shaping matrix as a continuous curve we can dynamically control, and use iterative dynamic programming along the nominal trajectory to find smoothly varying shaping matrices.

Contributions

In this work, we reformulate the reachable set computation problem as an optimal control problem. First, we propose a simple constraint-based set representation called a normotope, which is the sublevel set of a norm with a specified shaping matrix, center, and radius. We build a controlled embedding system for the original system using interval-based linear differential inclusions and logarithmic norms, with the observation that we may smoothly vary the shape matrix of the normotope as long as the offset term compensates accordingly—we call this the hypercontrol.

Next, we pose the following question: given an initial set for the system, can we find a hypercontrol mapping through the embedding space such that final reachable set has minimal volume? To address this problem, we define an optimal control problem to synthesize hypercontrol policies resulting in small terminal set volume. We show how a dynamic programming algorithms like iLQR can trade extra computations to iteratively optimize the accuracy of the reachable set estimate with respect to the hypercontrol. With two case studies, we show how (i) just a few cheap iterations can dramatically improve the reachable set accuracy, and (ii) that the algorithm can learn nontrivial policies which accumulate extra uncertainty earlier in the trajectory to shrink the final overapproximating reachable set.

II Background and Mathematical Preliminaries

II-A Notations

Let $\|\cdot\|$ denote a norm. Let $I$ denote the identity matrix of appropriate dimension and $\mathbf{1}$ denote the vector of all ones of appropriate dimension. For linear operator $A:\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{y}}$ between normed vector spaces, we set $\|A\|_{x\to y}=\sup_{x\in\mathbb{R}^{n_{x}}:\|x\|_{\mathbb{R}^{n_{x}}}=1}\|Ax\|_{\mathbb{R}^{n_{y}}}$ to be the operator norm of $A$. For linear operator $A:\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{x}}$, define the induced logarithmic norm $\mu(A)=\limsup_{h\searrow 0}\frac{\|I+hA\|_{x\to x}-1}{h}$. $GL(n)$ is the set of invertible $n\times n$ matrices. Let $\operatorname{co}\mathcal{S}$ denote the convex hull of $\mathcal{S}$.

II-B Nonlinear Systems and Reachable Sets

We consider nonlinear dynamical systems of the form

$$\dot{x}=f(t,x,w),$$ (1)

where $x\in\mathbb{R}^{n_{x}}$ is the state of the system, $w\in\mathbb{R}^{n_{w}}$ is a disturbance input to the system, and $f:\mathbb{R}_{\geq 0}\times\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{w}}\to\mathbb{R}^{n_{x}}$ is a $C^{1}$ smooth function.

Under these assumptions, the system (1) has a unique trajectory starting from any initial condition $x_{0}\in\mathbb{R}^{n_{x}}$ at $t_{0}$, under essentially bounded disturbance map $\mathbf{w}:[t_{0},\infty)\to\mathbb{R}^{n_{w}}$ defined for some neighborhood of $t_{0}$. Let $\phi_{f}(t;t_{0},x_{0},\mathbf{w})$ denote this trajectory.

Suppose we are given an initial set $\mathcal{X}_{0}\subseteq\mathbb{R}^{n_{x}}$ at initial time $t_{0}$, and a compact disturbance set $\mathcal{W}\subset\mathbb{R}^{n_{w}}$. Then we can define the (forward) reachable set of the system as

$$\mathcal{R}_{f}(t;t_{0},\mathcal{X}_{0},\mathcal{W})$$

$$:=\{\phi_{f}(t;t_{0},x_{0},\mathbf{w}):x_{0}\in\mathcal{X}_{0},\,\mathbf{w}:[t_{0},t_{f}]\to\mathcal{W}\}.$$

The reachable set is generally intractable to find in closed-form, both analytically and computationally, as it requires an infinite number of forward simulations of the system. The goal of this work and many others in the literature [20] is to efficiently compute a set overapproximation $\overline{\mathcal{R}}_{f}(t;t_{0},\mathcal{X}_{0},\mathcal{W})\supseteq\mathcal{R}_{f}(t;t_{0},\mathcal{X}_{0},\mathcal{W})$.

An overapproximation of the true reachable set can verify several safety specifications for the system—for instance, reach-avoid problems can be posed as

$$\underbrace{\overline{\mathcal{R}}(t_{f};t_{0},\mathcal{X}_{0},\mathcal{W})\subseteq\mathcal{G}}_{\text{reach $\mathcal{G}$ at $t=t_{f}$}},\ \underbrace{\overline{\mathcal{R}}(t;t_{0},\mathcal{X}_{0},\mathcal{W})\subseteq\mathcal{A}^{\complement}\ \forall t\in[t_{0},t_{f}]}_{\text{avoid $\mathcal{A}$ for every $t$}}.$$

In the coming sections, we propose a new set representation called a normotope, build a controlled embedding system whose trajectory provides a normotope overapproximating the true reachable set, and optimize over the choice of hypercontrol mappings in this embedding space.

III Normotopes: A Simple Constraint-Based Set Representation

In this section, we introduce the normotope, some basic properties, and a simple recipe to compute the normotope reachable set of a linear system using the adjoint dynamics.

III-A Definition and Examples

Definition 1.

Given a norm $\|\cdot\|$, a normotope is the set

$$\llbracket{\mathring{x}},\alpha,y\rrbracket:=\{x\in\mathbb{R}^{n_{x}}:\|\alpha(x-{\mathring{x}})\|\leq y\},$$

where ${\mathring{x}}\in\mathbb{R}^{n_{x}}$ is the center, $\alpha\in GL({n_{x}})$ is the square invertible ${n_{x}}\times{n_{x}}$ shape matrix, and $y>0$ is the offset.

Normotopes are a special case of a parametope, the general constraint-based set representation introduced in [12], and are dual to basic ellipsotopes [5] (letting $z=\tfrac{\alpha}{y}(x-{\mathring{x}})$, $\llbracket{\mathring{x}},\alpha,y\rrbracket=\{y\alpha^{-1}z+{\mathring{x}}:\|z\|\leq 1\}$, which is a basic ellipsotope with center ${\mathring{x}}$ and generator matrix $y\alpha^{-1}$).

Example 1 ($\ell_{2}$ normotopes).

The normotope

$$\llbracket{\mathring{x}},\alpha,y\rrbracket_{2} =\{x\in\mathbb{R}^{n_{x}}:\|\alpha(x-{\mathring{x}})\|_{2}\leq y\}$$

$$=\{x\in\mathbb{R}^{n_{x}}:(x-{\mathring{x}})^{T}\alpha^{T}\alpha(x-{\mathring{x}})\leq y^{2}\}$$

is equivalent to an ellipsoid centered at ${\mathring{x}}$, with weighting matrix $\alpha^{T}\alpha$ and radius $y^{2}$. $\alpha\in GL(n)$ implies $\alpha^{T}\alpha\succ 0$.

Example 2 ($\ell_{\infty}$ normotopes).

The normotope

$$\llbracket{\mathring{x}},\alpha,y\rrbracket_{\infty} =\{x\in\mathbb{R}^{n_{x}}:\|\alpha(x-{\mathring{x}})\|_{\infty}\leq y\}$$

$$=\{x\in\mathbb{R}^{n_{x}}:-y\mathbf{1}\leq\alpha(x-{\mathring{x}})\leq y\mathbf{1}\}$$

$$=\{x\in\mathbb{R}^{n_{x}}:\alpha{\mathring{x}}-y\mathbf{1}\leq\alpha x\leq\alpha{\mathring{x}}+y\mathbf{1}\}$$

is equivalent to the H-rep polytope given by $\{x:Hx\leq b\}$, $H=\begin{bmatrix}\alpha\\ -\alpha\end{bmatrix}$, $b=\begin{bmatrix}\alpha{\mathring{x}}+y\mathbf{1}\\ -\alpha{\mathring{x}}+y\mathbf{1}\end{bmatrix}$, with $2n$ faces and $2^{n}$ vertices.

Example 3 ($\ell_{1}$ normotopes).

The normotope

$$\llbracket{\mathring{x}},\alpha,y\rrbracket_{1} =\{x\in\mathbb{R}^{n_{x}}:\|\alpha(x-{\mathring{x}})\|_{1}\leq y\}$$

$$=\{x\in\mathbb{R}^{n_{x}}:-y\mathbf{1}\leq S\alpha(x-{\mathring{x}})\leq y\mathbf{1}\}$$

$$=\{x\in\mathbb{R}^{n_{x}}:S\alpha{\mathring{x}}-y\mathbf{1}\leq S\alpha x\leq S\alpha{\mathring{x}}+y\mathbf{1}\}$$

where $S\in\mathbb{R}^{2^{n_{x}}\times{n_{x}}}$ is the matrix consisting of all $2^{n_{x}}$ possible choices of signs, i.e., the $k$-th row $S_{k}=(s^{k}_{1},\dots,s^{k}_{n})$ where $s^{k}_{j}\in\{+1,-1\}$. is equivalent to the H-rep polytope given by $\{x:Hx\leq b\}$, $H=\begin{bmatrix}S\alpha\\ -S\alpha\end{bmatrix}$, $b=\begin{bmatrix}S\alpha{\mathring{x}}+y\mathbf{1}\\ -S\alpha{\mathring{x}}+y\mathbf{1}\end{bmatrix}$, with $2^{n}$ faces and $2n$ vertices.

Remark 1 (Nomenclature).

We choose the name “normotope” over “norm-ball” to emphasize how we expect the shape matrix $\alpha$ to vary along reachable set trajectories.

III-B The LTV Case

In the undisturbed linear time-varying case (LTV) with a normotope initial set, the adjoint equation provides the true reachable set of the system as a normotope with fixed offset.

Proposition 1.

For the LTV system

$$\dot{x}(t)=A(t)x(t),$$

the true reachable set is given by

$$\mathcal{R}_{f}(t;t_{0},\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket)=\llbracket{\mathring{x}}(t),\alpha(t),y_{0}\rrbracket,$$

where $t\mapsto{\mathring{x}}(t),\alpha(t)$ satisfy the following ODEs,

$$\dot{{\mathring{x}}}(t) =A(t){\mathring{x}}(t),\quad\dot{\alpha}(t)=-\alpha(t)A(t).$$

with initial conditions ${\mathring{x}}(t_{0})={\mathring{x}}_{0}$, $\alpha(t_{0})=\alpha_{0}$.

Proof.

We first recall how the flow map is given by $\phi_{f}(t;t_{0},x_{0})=\Phi(t;t_{0})x_{0}$, where $\Phi$ is the state transition matrix satisfying the matrix ODE

$$\dot{\Phi}(t;t_{0})=A(t)\Phi(t;t_{0}),\quad\Phi(t_{0};t_{0})=\mathbf{I},$$

with the property $\Phi^{-1}(t;t_{0})=\Phi(t_{0};t)$. Thus, the true reachable set is given by

$$\mathcal{R}(t;t_{0},\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket)=\{\Phi(t;t_{0})x:\|\alpha_{0}(x-{\mathring{x}}_{0})\|\leq y_{0}\}$$

$$=\{z:\|\alpha_{0}(\Phi(t_{0};t)z-{\mathring{x}}_{0})\|\leq y_{0}\}$$

$$=\{z:\|\alpha_{0}\Phi(t_{0};t)(z-\Phi(t;t_{0}){\mathring{x}}_{0})\|\leq y_{0}\}$$

$$=\llbracket\Phi(t;t_{0}){\mathring{x}},\alpha_{0}\Phi(t_{0};t),y_{0}\rrbracket=\llbracket{\mathring{x}}(t),\alpha(t),y_{0}\rrbracket,$$

where $\alpha(t)$ satisfies the following,

$$\dot{\alpha}(t)=\alpha_{0}\frac{\mathrm{d}}{\mathrm{d}t}\Phi(t;t_{0})^{-1}=-\alpha_{0}\Phi(t;t_{0})^{-1}\dot{\Phi}(t;t_{0})\Phi(t;t_{0})^{-1}$$

$$=-\alpha_{0}\Phi(t_{0};t)A(t)\Phi(t;t_{0})\Phi(t;t_{0})^{-1}=-\alpha(t)A(t),$$

and ${\mathring{x}}(t)$ satisfies $\dot{{\mathring{x}}}(t)=A(t){\mathring{x}}(t)$. ∎

We can interpret the ODE $\dot{\alpha}(t)=-\alpha(t)A(t)$ as each row $\alpha_{j}^{T}$ evolving according to the adjoint dynamics $\dot{\alpha}_{j}(t)=-A(t)^{T}\alpha_{j}(t)$.

In the nonlinear systems case, the true reachable set is generally not itself a normotope. In the next section, we demonstrate how a similar embedding system approach can obtain guaranteed overapproximating normotope reachable sets by adding an appropriate dynamic compensation term to the offset $y$ using logarithmic and induced matrix norms.

IV Normotope Embeddings

In Proposition 1, we showed that the true reachable set in the linear time-varying case was given by a single forward simulation of another ODE (what we will refer to as an embedding system), by evolving the center ${\mathring{x}}$ using the system dynamics and the rows of the shape matrix $\alpha$ using the adjoint dynamics. In this section, for nonlinear systems, we demonstrate how a similar controlled embedding approach can be used to bound the behavior of any nonlinear dynamics.

IV-A Dynamics Abstraction: Linear Differential Inclusions

A linear differential inclusion (LDI) [21, p.51] encompassing the error dynamics of the system (1) is an inclusion where for every $x\in\mathcal{X}$ and $w\in\mathcal{W}$, for prespecified ${\mathring{x}}\in\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$ and ${\mathring{w}}\in\mathcal{W}\subseteq\mathbb{R}^{n_{w}}$,

$$f(t,x,w)-f(t,{\mathring{x}},{\mathring{w}})\in\mathcal{M}^{x}(t)(x-{\mathring{x}})+\mathcal{M}^{w}(t)(w-{\mathring{w}}),$$ (2)

where $\mathcal{M}^{x}(t)\subseteq\mathbb{R}^{{n_{x}}\times{n_{x}}}$ and $\mathcal{M}^{w}(t)\subseteq\mathbb{R}^{{n_{x}}\times{n_{w}}}$ are sets of matrices. There are many ways to obtain matrix sets satisfying (2). For example, if $\overline{\operatorname{co}}\{\frac{\partial f}{\partial x}(t,x,w)\}\subseteq\mathcal{M}^{x}(t)$ and $\overline{\operatorname{co}}\{\frac{\partial f}{\partial w}(t,x,w)\}\subseteq\mathcal{M}^{w}(t)$ for every $(x,w)\in\mathcal{X}\times\mathcal{W}$ (where $\overline{\operatorname{co}}$ is the closed convex hull), then (2) holds for every $(x,w)\in\mathcal{X}\times\mathcal{W}$ [17, Prop. 1], [21, p.55].

We briefly discuss an algorithmic method from [17, Cor. 1] to obtain the matrices $\mathcal{M}^{x}$ and $\mathcal{M}^{w}$ using interval analysis on the first partial derivatives of $f$. Suppose we are bounding a $C^{1}$ function $g:\mathbb{R}^{n}\to\mathbb{R}^{m}$, and we are given given an interval set $Z_{1}\times\cdots\times Z_{n}\subseteq\mathbb{R}^{n}$. Then the following implication holds, for fixed $\mathring{z}\in Z_{1}\times\cdots\times Z_{n}$,

$$\frac{\partial g_{i}}{\partial x_{j}}(Z_{1},\dots,Z_{j},\mathring{z}_{j+1},\dots,\mathring{z}_{n})\subseteq[\mathcal{M}]_{ij}$$

$$\implies g(z)-g(\mathring{z})\in[\mathcal{M}](z-\mathring{z}),$$

for every $z\in Z_{1}\times\cdots Z_{n}$ [17, Cor. 1], where the RHS is evaluated using interval matrix multiplication [22]. We can apply this procedure to the map $f(t,\cdot,\cdot):\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{w}}\to\mathbb{R}^{n_{x}}$ by defining an equivalent map $\hat{f}_{t}:\mathbb{R}^{{n_{x}}+{n_{w}}}\to\mathbb{R}^{n_{x}}$ to obtain the bound

$$f(t,x,w)-f(t,{\mathring{x}},{\mathring{w}})\in[\mathcal{M}(t)]\begin{bmatrix}x-{\mathring{x}}\\ w-{\mathring{x}}\end{bmatrix},$$

and finally, we can split $[\mathcal{M}(t)]\subseteq\mathbb{R}^{{n_{x}}\times({n_{x}}+{n_{w}})}$ into two interval matrices $[\mathcal{M}^{x}(t)]\subseteq\mathbb{R}^{{n_{x}}\times{n_{x}}}$ and $[\mathcal{M}^{w}(t)]\subseteq\mathbb{R}^{{n_{x}}\times{n_{w}}}$.

We simplify the expression by replacing the interval matrices $[\mathcal{M}^{x}(t)]$ and $[\mathcal{M}^{w}(t)]$ with the convex hull of a finite set of corners. For example, if the interval matrix $[\mathcal{M}^{x}(t)]$ has $4$ non-singleton entries, then there are $2^{4}=16$ corners $M_{i}^{x}(t)$ to consider. Thus, we obtain a bound of the following form,

$$f(t,x,w)- f(t,{\mathring{x}},{\mathring{w}})\in$$

$$\operatorname{co}\{M_{i}^{x}(t)\}_{i}(x-{\mathring{x}})+\operatorname{co}\{M_{j}^{w}(t)\}_{j}(w-{\mathring{w}}).$$

As we will see in Theorem 1 below, these finite sets will help us build a simple embedding system by evaluating logarithmic and induced matrix norms over each corner.

IV-B Controlled Embedding System

In Theorem 1, we apply build a particularly well structured embedding system for the special case of a normotope set with a linear differential inclusion bounding the error dynamics to the center of the normotope.

Theorem 1 (Controlled normotope embedding).

Let $\|\cdot\|$ be a differentiable norm, the $\ell_{1}$-norm, or the $\ell_{\infty}$-norm. Let $\mu$ be the induced logarithmic norm. Define the following controlled embedding system,

$$\dot{{\mathring{x}}}(t) =f(t,{\mathring{x}},{\mathring{w}}),$$ (3a)

$$\dot{\alpha}(t) =U(t,{\mathring{x}},\alpha,y),$$ (3b)

$$\begin{split}\dot{y}(t)&=\left(\max_{i}\mu(U(t,{\mathring{x}},\alpha,y)\alpha^{-1}+\alpha M^{x}_{i}(t)\alpha^{-1})\right)y\\ &\quad\quad+\max_{j}\|\alpha M^{w}_{j}(t)\|_{w\to x},\end{split}$$ (3c)

where $f(t,x,w)-f(t,{\mathring{x}},{\mathring{w}})\in\operatorname{co}\{M^{x}_{i}(t)\}(x-{\mathring{x}})+\operatorname{co}\{M^{w}_{i}(t)\}(w-{\mathring{w}})$ for every $(x,w)\in\partial\llbracket{\mathring{x}},\alpha,y\rrbracket\times\mathcal{W}$. Then for every $t\geq 0$,

$$\mathcal{R}_{f}(t,\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket,\mathcal{W})\subseteq\llbracket{\mathring{x}}(t),\alpha(t),y(t)\rrbracket,$$

where $t\mapsto({\mathring{x}}(t),\alpha(t),y(t))$ is the trajectory of the system (3) under any locally Lipschitz map $U$.

Proof.

((i) $\|\cdot\|$ is differentiable). For $(x,w)\in\partial\llbracket{\mathring{x}},\alpha,y\rrbracket\times\mathcal{W}$, there is $M^{x}\in\operatorname{co}\{M^{x}_{i}\}_{i}$ and $M^{w}\in\operatorname{co}\{M^{w}_{j}\}_{j}$ such that

$$\xi(t,x,w):=\partial_{\alpha}g[U]+\partial_{x}g(M^{x}(x-{\mathring{x}})+M^{w}(w-{\mathring{w}}))$$

$$=\tfrac{\partial\|\cdot\|}{\partial z}\big|_{\alpha(x-{\mathring{x}})}(U(x-{\mathring{x}})+\alpha(M^{x}(x-{\mathring{x}})+M^{w}(w-{\mathring{w}})))$$

$$=\lim_{h\to 0}\tfrac{\|\alpha(x-{\mathring{x}})+hU(x-{\mathring{x}})+hM^{x}(x-{\mathring{x}})+hM^{w}(w-{\mathring{w}})\|-\|\alpha(x-{\mathring{x}})\|}{h}$$

$$\leq\limsup_{h\searrow 0}\tfrac{\|(I+h(U\alpha^{-1}+\alpha M^{x}\alpha^{-1}))\alpha(x-{\mathring{x}})\|-\|\alpha(x-{\mathring{x}})\|}{h}$$

$$\quad+\|\alpha M^{w}(w-{\mathring{w}})\|$$

$$\leq\limsup_{h\searrow 0}\tfrac{\|I+h(U\alpha^{-1}+\alpha M^{x}\alpha^{-1})\|_{x\to x}-1}{h}\|\alpha(x-{\mathring{x}})\|$$

$$\quad+\|\alpha M^{w}(w-{\mathring{w}})\|$$

$$=\mu(U\alpha^{-1}+\alpha M^{x}\alpha^{-1})y+\|\alpha M^{w}\|_{\operatorname{op}}$$

$$\leq\max_{i}\mu(U\alpha^{-1}+\alpha M^{x}_{i}\alpha^{-1})y+\max_{j}\|\alpha M^{w}_{j}\|_{w\to x}.$$

where the final inequality holds since $\mu$ and $\|\cdot\|_{w\to x}$ are convex functions [23, Lemma 2.3] with arguments linear in $M^{x}$ and $M^{w}$ respectively. [12, Cor. 1] concludes the proof.

((ii) $\|\cdot\|=\|\cdot\|_{\infty}$) We first note that the normotope $\{\|\alpha(x-{\mathring{x}})\|_{\infty}\leq y\}$ is equivalent to the symmetric H-rep polytope $\{-\hat{y}\leq\alpha(x-{\mathring{x}})\leq\hat{y}\}$ where $\hat{y}=y\mathbf{1}\in\mathbb{R}^{n_{x}}$. Then for any $x$ satisfying $-\hat{y}\leq\alpha(x-{\mathring{x}})\leq\hat{y}$ and $|\alpha_{k}(x-{\mathring{x}})|=y$ for some $k$ (where $\alpha_{k}$ is the $k$-th row of $\alpha$),

$$\xi_{k}(t,x,w):=\dot{\alpha}_{k}(x-{\mathring{x}})+\alpha_{k}(f(t,x,w)-f(t,{\mathring{x}},{\mathring{w}}))$$

$$=\dot{\alpha}_{k}(x-{\mathring{x}})+\alpha_{k}(M^{x}(x-{\mathring{x}})+M^{w}(w-{\mathring{w}}))$$

$$=(\dot{\alpha}_{k}\alpha^{-1}+\alpha_{k}M^{x}\alpha^{-1})\alpha(x-{\mathring{x}})+\alpha_{k}M^{w}(w-{\mathring{w}})$$

$$\leq\big((\dot{\alpha}_{k}\alpha^{-1}+\alpha_{k}M^{x}\alpha^{-1})_{k}$$

$$\quad+\textstyle\sum_{j\neq k}|(\dot{\alpha}_{k}\alpha^{-1}+\alpha_{k}M^{x}\alpha^{-1})_{j}|\big)y+\alpha_{k}M^{w}(w-{\mathring{w}})$$

$$\leq\mu_{\infty}(U\alpha^{-1}+\alpha M^{x}\alpha^{-1})y+\|\alpha M^{w}\|_{w\to x}$$

$$\leq\max_{i}\mu_{\infty}(U\alpha^{-1}+\alpha M^{x}_{i}\alpha^{-1})y+\max_{j}\|\alpha M^{w}_{j}\|_{w\to x}$$

by definition of the $\mu_{\infty}$ log norm [23, Table 2.1], and concluding with the same logic as part (i). The resulting H-polytope embedding recovers the same trajectory as (3) with the equivalence $\hat{y}=y\mathbf{1}$.

((iii) $\|\cdot\|=\|\cdot\|_{1}$) We first note that

$$\|x\|_{1} =\textstyle\sum_{i}|x_{i}|=\max_{s_{j}\in S}(s_{j}^{T}x)=\|Sx\|_{\infty},$$

where $S\in\mathbb{R}^{2^{n_{x}}\times{n_{x}}}$ is the matrix consisting of all $2^{n_{x}}$ possible choices of signs, i.e., the $k$-th row $S_{k}=(s^{k}_{1},\dots,s^{k}_{n})$ where $s^{k}_{j}\in\{+1,-1\}$. Thus, the normotope $\{\|\alpha(x-{\mathring{x}})\|_{1}\leq y\}$ is equivalent to the symmetric H-rep polytope given by $\{x:-\hat{y}\leq S\alpha(x-{\mathring{x}})\leq\hat{y}\}$, where $\hat{y}=y\mathbf{1}$. Let $S^{+}$ be a left inverse of $S$ ($S^{+}S=\mathbf{I}$). Then for any $x$ satisfying $-\hat{y}\leq S\alpha(x-{\mathring{x}})\leq\hat{y}$ and $|S_{k}\alpha(x-{\mathring{x}})|=y$ for some $k$ (where $\alpha_{k}$ is the $k$-th row of $\alpha$),

$$\xi_{k}(t,x,w)=(S_{k}\dot{\alpha})(x-{\mathring{x}})$$

$$\quad+S_{k}\alpha(M^{x}(x-{\mathring{x}})+M^{w}(w-{\mathring{w}}))$$

$$=(S_{k}\dot{\alpha}\alpha^{-1}S^{+}+S_{k}\alpha M^{x}\alpha^{-1}S^{+})S\alpha(x-{\mathring{x}})$$

$$\quad+S^{k}\alpha M^{w}(w-{\mathring{w}})$$

$$\leq\mu_{\infty}(S(U\alpha^{-1}+\alpha M^{x}\alpha^{-1})S^{+})y+\|S\alpha M^{w}\|_{\infty\to w},$$

concluding with the logic of part (ii). Finally, for any $A$,

$$\|SAS^{+}\|_{\infty} =\sup_{z:\|z\|_{\infty}=1}\|SAS^{+}z\|_{\infty}=\sup_{x:\|Sx\|_{\infty}=1}\|SAx\|_{\infty}$$

$$=\sup_{x:\|x\|_{1}}\|Ax\|_{1}=\|A\|_{1},$$

thus, plugging into their definitions, we see that $\mu_{\infty}(SAS^{+})=\mu_{1}(A)$ and $\|SA\|_{w\to\infty}=\|A\|_{w\to 1}$. ∎

In practice, we first overapproximate the normotope $\llbracket{\mathring{x}}(t),\alpha(t),y(t)\rrbracket$ with an interval subset, then apply the methodology from the previous section to obtain the matrices $\{M_{i}^{x}(t)\}_{i}$ and $\{M_{j}^{w}(t)\}_{j}$. Finally, we note that the dynamics (3c) are not the only possible choice of embedding dynamics; in fact, any expression satisfying the hypotheses of [12, Cor. 1] would constitute a valid embedding. However, when the dynamics are abstracted using a linear differential inclusion as (2), the expression (3c) is a natural choice with direct connections to contraction analysis [24, 19, 23].

Remark 2 (Connection to contraction theory).

The $U=0$ case of Theorem 1 recovers the fixed norm results from [13, 16, 17] for the norm with fixed shape matrix $\|x\|_{\alpha}=\|\alpha x\|$. The introduction of the hypercontrol $U$ allows us to continuously vary the shape matrix along the trajectory of the normotope embedding system.

V An Optimal Control Perspective Towards Reachable Set Computation

In this section, we synthesize the policy $U$ from (3) using an optimal control formulation. For clarity of presentation, we omit the disturbance input in this section and consider the nonlinear system

$$\dot{x}=f(t,x),$$

where $f:\mathbb{R}\times\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{x}}$ is $C^{1}$ differentiable, but we note that all results in this section can easily be generalized to the case with a compact disturbance.

As seen in Proposition 1, for LTV systems the best choice of embedding dynamics is $\dot{\alpha}=-\alpha A(t)$. We choose hypercontrol policies of the form

$$U(t,{\mathring{x}},\alpha,y)=-\alpha\frac{\partial f}{\partial x}({\mathring{x}})+U_{\text{ff}}(t),$$

where $U_{\text{ff}}$ is a feed-forward input to the embedding system. As discussed in [25], closing the loop with the adjoint of the linearization cancels the first order expansion of the dynamics and generally provides reasonable results when the offset is small and the system is close to linear. In the rest of this section, using the controlled embedding system from the previous section, we formulate an optimal control problem and use iterative dynamic programming to optimize the accuracy of our reachable set overapproximations.

V-A Cost Function Design

Suppose we are given an initial set $\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket$ at $t=t_{0}$, and we are interested in verifying a goal specification at $t=t_{f}$. For a goal specification, we are purely interested in the size of the reachable set at the final time, and therefore, we design an optimal control problem with a pure terminal cost. The following Lemma recalls a standard result whereby maximizing the log-determinant of the shape matrix minimizes the volume of the parametope set.

Lemma 1.

For any norm $\|\cdot\|$, a minimizer of the objecive function

$$\Phi({\mathring{x}},\alpha,y)=-\log\det(\alpha^{T}\alpha/y^{2})$$ (4)

minimizes the volume of the normotope set

$$\operatorname{vol}(\llbracket{\mathring{x}},\alpha,y\rrbracket).$$

Proof.

The normotope set

$$\llbracket{\mathring{x}},\alpha,y\rrbracket =\{x:\|\alpha(x-{\mathring{x}})\|\leq y\}=\{x:\|\tfrac{\alpha}{y}(x-{\mathring{x}})\|\leq 1\}$$

$$=\{(\tfrac{\alpha}{y})^{-1}z+{\mathring{x}}:\|z\|\leq 1\}$$

is an affine image of the norm ball $\mathcal{B}_{1}(0)$. Thus, the volume

$$\operatorname{vol}(\llbracket{\mathring{x}},\alpha,y\rrbracket)=|\det((\tfrac{\alpha}{y})^{-1})|\operatorname{vol}(\mathcal{B}_{1}(0))$$

is inversely proportional to $|\det(\tfrac{\alpha}{y})|$. Since $\log((\cdot)^{2})$ is a monotone function for positive inputs, we can maximize

$$\log(|\det(\tfrac{\alpha}{y})|^{2})=\log(\det(\tfrac{\alpha^{T}}{y})\det(\tfrac{\alpha}{y}))=\log\det(\tfrac{\alpha^{T}\alpha}{y^{2}}),$$

completing the proof. ∎

We implement normotopes and the terminal cost using JAX for Python, taking care to respect the stringent requirements of JAX, which provides us several crucial capabilities. Normotopes are natively implemented as a Pytree node class, allowing us to use them as a valid JAX type for operations such as automatic differentiation. For the terminal cost (4), we use jnp.linalg.cholesky instead of composing log and det for numerical stability.

V-B Open-Loop Hypercontrol Synthesis: Iterative Linear Quadratic Regulator

Algorithm 1 Reach-iLQR

Input: initial set $X_{0}=\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket$, time interval $[t_{0},t_{f}]$
Parameters: step size $\gamma>0$, regularizer $R\succ 0$, volume threshold $\Phi_{\text{max}}\in\mathbb{R}$

1: $U_{\text{ff}}^{(0)}(t)\leftarrow 0$, $d(t)\leftarrow 0$, $K(t)\leftarrow 0$

2: $i\leftarrow 1$

3: repeat

4:  Forward simulate under the updated control $U_{\text{ff}}^{(i)}$ from (9) to obtain $X^{(i)}(t)$ until $\Phi>\Phi_{\text{max}}$ or $t=t_{f}$.

5:  Compute terminal conditions (7) using automatic differentiation on terminal cost $\Phi$ from (4).

6:  Backward simulate Ricatti equations (6) using automatic differentiation to obtain $d(t),K(t)$ from (8).

7:  $i\leftarrow i+1$

8: until termination criterion reached

Suppose we are given an initial set $\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket$. In this section, we apply the continuous-time iterative linear quadratic regulator approach (iLQR) [26] to numerically optimize the following optimal control problem.

$$\operatornamewithlimits{minimize}_{U_{\text{ff}}:[t_{0},t_{f}]\to\mathbb{R}^{{n_{x}}\times{n_{x}}}}\ \Phi({\mathring{x}}(t_{f}),\alpha(t_{f}),y(t_{f}))$$ (5)

$$\text{s.t. }\quad f(t,x)-f(t,{\mathring{x}})\in\operatorname{co}\{M_{i}^{x}(t)\}(x-{\mathring{x}}),$$

$$\begin{aligned} \dot{{\mathring{x}}}(t)&=f(t,{\mathring{x}}),\\ \dot{\alpha}(t)&=-\alpha\tfrac{\partial f}{\partial x}(t,{\mathring{x}})+U_{\text{ff}}(t)\\ \dot{y}(t)&=\left(\max_{i}\mu(U_{\text{ff}}(t)\alpha^{-1}+\alpha(M_{i}^{x}(t)-\tfrac{\partial f}{\partial x}(t,{\mathring{x}}))\alpha^{-1})\right)y\\ \end{aligned}$$

$$({\mathring{x}}(t_{0}),\alpha(t_{0}),y(t_{0}))=({\mathring{x}}_{0},\alpha_{0},y_{0})$$

For the rest of this section, we write $X=\operatorname{vec}({\mathring{x}}_{0},\alpha_{0},y_{0})\in\mathbb{R}^{n_{X}}=\mathbb{R}^{{n_{x}}+n_{x}^{2}+1}$, and we abbreviate the embedding dynamics (3) as $\dot{X}=F(t,X,U)$.

Next, we describe Algorithm 1, which uses continuous-time iterative linear quadratic regulation (iLQR) to improve the feedforward term $U_{\text{ff}}$. Suppose we are at timestep $i$, given $U_{\text{ff}}^{(i)}(t)$ and the corresponding embedding system trajectory $X^{(i)}(t)$ from forward simulating the embedding dynamics (3).

V-B1 Backward Pass

Since we have no running cost, the backwards pass simplifies to the following ODEs,

$$-\dot{\sigma}(t) =-\tfrac{1}{2}Q_{U}^{T}(t)R^{-1}Q_{U}(t),$$ (6a)

$$-\dot{s}(t) =F_{X}^{T}(t)s(t)-Q_{UX}^{T}(t)R^{-1}Q_{U}(t),$$ (6b)

$$-\dot{S}(t) =Q_{XX}(t)-Q_{UX}^{T}(t)R^{-1}Q_{UX}(t),$$ (6c)

where $R\succ 0$ is some regularization matrix,

$$Q_{U}(t) =F_{U}^{T}(t)s(t),$$

$$Q_{XX}(t) =F_{X}^{T}(t)S(t)+S(t)F_{X}(t),$$

$$Q_{UX}(t) =F_{U}^{T}(t)S(t),$$

and the following

$$F_{U}^{T}(t) =\frac{\partial F}{\partial U}(t,X^{(i)}(t),U_{\text{ff}}^{(i)}(t)),$$

$$F_{X}^{T}(t) =\frac{\partial F}{\partial X}(t,X^{(i)}(t),U_{\text{ff}}^{(i)}(t)),$$

are computed using automatic differentiation in JAX. The ODEs are initialized by automatically differentiating the terminal cost $\Phi$ with respect to the state $X$ as

$$\sigma(t_{f}) =\Phi(X(t_{f})),$$ (7a)

$$s(t_{f}) =\frac{\partial\Phi}{\partial X}(X(t_{f}))^{T},$$ (7b)

$$S(t_{f}) =\frac{\partial^{2}\Phi}{\partial X^{2}}(X(t_{f})).$$ (7c)

Finally, we obtain curves $d(t)$ and $K(t)$ as

$$d(t)=R^{-1}Q_{U}(t),\quad K(t)=R^{-1}Q_{UX}(t).$$ (8)

The key result from iLQR is that the variational system (first order expansion of the dynamics)

$$\dot{\delta X}=F_{X}^{T}(t)\delta X(t)+F_{U}^{T}(t)\delta U(t)$$

and the second order expansion of the value function

$$V(x,t)\approx\sigma(t)+s(t)^{T}\delta X(t)+\tfrac{1}{2}\delta X(t)^{T}S(t)\delta X(t)$$

satisfy the Hamilton-Jacobi-Bellman PDE with feedback controller $\delta U(t)=-d(t)-K(t)\delta X(t)$ (for the full derivation, please see [26]).

V-B2 Forward Pass

From the backward pass, we obtain the curves $d(t)$ and $K(t)$. We then perform another forward sweep by updating $U_{\text{ff}}$ using the following (where $i$ has been incremented),

$$U_{\text{ff}}^{(i)}(t)=U_{\text{ff}}^{(i-1)}(t)-\gamma d(t)-K(t)(X^{(i)}(t)-X^{(i-1)}(t)),$$ (9)

where $\gamma>0$ is the step size, $X^{(i-1)}$ is the previous embedding trajectory under $U_{\text{ff}}^{(i-1)}$, and $X^{(i)}$ is the new embedding trajectory under $U_{\text{ff}}^{(i)}$. Theoretically, when $\gamma$ is small enough, the displacement $\delta X=X^{(i)}-X^{(i-1)}$ along the trajectory is small enough where the first order dynamics expansion and second order value expansion remain accurate.

While the forward pass always produces a valid reachable set overapproximation, it is often possible that the reachable set blows up before the prescribed final time $t_{f}$. To combat this situation, we introduce a threshold $\Phi_{\text{max}}$, where we stop the forward simulation if $\Phi(X(t))>\Phi_{\text{max}}$. When this happens, we proceed as if the final time $t_{f}$ was the time instance before the threshold was violated; in practice, we find that iterates will eventually tend towards policies which reach the final time. The result is an algorithm where each iteration provides a valid overapproximate reachable tube, tending towards those with smaller terminal set volumes.

Remark 3 (Comparison to the literature).

Previous approaches [13, 16, 17] break the time interval $[t_{0},t_{f}]$ into segments, solve intermediate semidefinite programs to optimize over the potential choices of norms along each individual segment, and ensure the reachable set containment by overapproximating the norm ball at the end of the previous segment with the new norm ball at the start of the next segment. Algorithm 1 instead provides a continuous formulation where the hypercontrol smoothly varies the shaping matrix along the nominal trajectory, both avoiding successive norm overapproximations and allowing us to backpropogate the terminal set volume all the way to the initial timestep.

VI Examples

In the following examples, we demonstrate how Algorithm 1 can (i) be used to quickly modify the reachable set to obtain drastically improve estimates in a short amount of time, and (ii) learn nontrivial open-loop policies for more difficult benchmarks. For all ODE simulations, we use Euler integration with a step size of $0.01$.

Figure 1: Top: Overapproximating reachable sets for the robot arm for $t\in[0,10]$ are plotted for various iteration numbers of Algorithm 1, with Monte Carlo samples in red. Bottom: The terminal cost $\Phi$ is plotted as a function of the iLQR iteration number in blue, while the cumulative runtime is in red. After 5 iterations of iLQR, the overapproximate reachable set reaches the final timestep. Iterating further, we can trade runtime for less overconservatism in the final reachable set.

Figure 2: Top: Overapproximating reachable sets for the Van der Pol oscillator for $t\in[0,7]$ are shown for various iteration numbers of Algorithm 1, with Monte Carlo samples in red. Bottom: The function $\Phi(X(t))$ is plotted as a function of time for the same iteration numbers. While the pure adjoint dynamics ($i=1$) does well for small horizons, the error quickly accumulates over longer horizons and cannot reach the goal of $t_{f}=7$s. On the other hand, after 3000 iterations of iLQR, the algorithm qualitatively learns to accumulate extra error in earlier timesteps so it can reach and minimize the final timestep.

VI-A Robot Arm

Consider the dynamics of a robot arm from [27, 15]

$$\dot{q}_{1} =z_{1},\quad\dot{q}_{2}=z_{2},$$

$$\dot{z}_{1} =\tfrac{1}{mq_{2}^{2}+ML^{2}/3}(-2mq_{2}z_{1}z_{2}-k_{d_{1}}z_{1}+k_{p_{1}}(u_{1}-q_{1})),$$

$$\dot{z}_{2} =q_{2}z_{1}^{2}+\tfrac{1}{m}(-k_{d_{2}}z_{2}+k_{p_{2}}(u_{2}-q_{2})),$$

with $u_{1}=2$, $u_{2}=1$, $m=M=1$, $L=\sqrt{3}$, $k_{p_{1}}=2$, $k_{p_{2}}=1$, $k_{d_{1}}=2$, $k_{d_{2}}=1$. We use the initial set $\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket_{2}$, where $\alpha_{0}=P^{1/2}/0.1$ with $P$ and ${\mathring{x}}_{0}$ from [15], and $y_{0}=1$. We note that this initial set has a radius $2.5$ times larger than the largest radius considered in the comparison example from [17], and thus a volume that is $39$ times larger. We run Algorithm 1 with $t_{f}=10$, $\gamma=1$, $R=20$, and $\Phi_{\text{max}}=-0.1$ for 20 iterations, taking a total of 7.5 seconds, and the results are shown in Figure 1.

VI-B Van der Pol Oscillator

Consider the Van der Pol oscillator written in usual state space coordinates as

$$\dot{x}_{1}=x_{2},\quad\dot{x}_{2}=-x_{1}+\eta(1-x_{1}^{2})x_{2},$$

with nonlinearity $\eta=1$. We consider the initial set $\llbracket{\mathring{x}}_{0},\alpha_{0},y_{0}\rrbracket_{2}$ for ${\mathring{x}}_{0}=[-2\ 0]^{T}$, $\alpha_{0}=80I$, and $y_{0}=1$.

We split the iLQR iterations into two phases. We first start with a small regularization of $R=0.5I$ as an exploratory phase, keeping track of the best reachable set computation in all the iterates (the one that reaches the furthest before violation $\Phi>\Phi_{\text{max}}=-1.75$, with the smallest volume). After 750 iterations, we revert to the best iterate from the first phase ($i=707$) and use a larger regularizer $R=5I$ for another 750 iterations to help reach the end. All 1500 iterations took a total of 23.8 seconds to execute.

VII Conclusion

In this work, we presented an optimal control framework for synthesizing guaranteed overapproximate reachable sets with small terminal volume using iterative dynamic programming. Our framework builds off the controlled parametric embedding approach described in [12], constructing a special case for normotope sets using logarithmic and induced matrix norms. In future work, we plan to use techniques from reinforcement learning to learn closed-loop policies using offline training to mimic the output of Algorithm 1, which would be quick to evaluate for online reachable set computations.

Appendix

VII-A The LTV Case: Adjoint Satisfies Pontryagin’s Principle

As a sanity check, we verify that the adjoint equation satisfies Pontryagin’s Minimum Principle (PMP) for the problem (5) in the LTV case. Consider the linear system

$$\dot{x}=A(t)x(t),$$

and the embedding system from Theorem 1, simplifying to

$$\dot{{\mathring{x}}}=A(t){\mathring{x}}(t),\ \dot{\alpha}=-\alpha(t)A(t)+\tilde{U}(t)\alpha(t),\ \dot{y}=\mu(\tilde{U}(t))y,$$

where we have made the variable substitution $\tilde{U}(t)=(\alpha(t)A(t)+U(t))\alpha(t)^{-1}$. Consider the volume terminal cost from (4),

$$\Phi({\mathring{x}},\alpha,y)=-\log\det(\alpha^{T}\alpha)+2n\log y.$$

The Hamiltonian of this problem is

$$H(t,{\mathring{x}},\alpha,y,p_{\mathring{x}},p_{\alpha},p_{y},\tilde{U})$$

$$=\langle p_{\mathring{x}},A{\mathring{x}}\rangle-\langle p_{\alpha},\alpha A\rangle+\langle p_{\alpha},\tilde{U}\alpha\rangle+\tfrac{1}{2}p_{y}y\lambda_{\max}(\tilde{U}^{T}+\tilde{U}),$$

thus, the costate dynamics satisfy

$$\dot{p}_{\alpha} =-H_{\alpha}=p_{\alpha}A^{T}-\tilde{U}^{T}p_{\alpha}$$

$$\dot{p}_{y} =-H_{y}=-\tfrac{1}{2}p_{y}\lambda_{\max}(\tilde{U}^{T}+\tilde{U})$$

with terminal costates given by

$$p_{\alpha}(t_{f})=\tfrac{\partial\Phi}{\partial\alpha}=-2\alpha(t_{f})^{-T},\quad p_{y}(t_{f})=\tfrac{\partial\Phi}{\partial y}=2n/y(t_{f}).$$

Set $\tilde{U}^{*}=0$. Then $\dot{y}=0$, $\dot{p}_{y}=0$, and

$$\frac{\mathrm{d}}{\mathrm{d}t}\alpha p_{\alpha}^{T}=\dot{\alpha}p_{\alpha}^{T}+\alpha\dot{p}_{\alpha}^{T}=-\alpha Ap_{\alpha}^{T}+\alpha Ap_{\alpha}^{T}=0$$

$$\implies\alpha(t)p_{\alpha}(t)^{T}=\alpha(t_{f})p_{\alpha}(t_{f})^{T}=-2\mathbf{I}$$

Plugging into the Hamiltonian,

$$H(\dots,\tilde{U})=\langle p_{\mathring{x}},A{\mathring{x}}\rangle-\langle p_{\alpha},\alpha A\rangle+\langle p_{\alpha},\tilde{U}\alpha\rangle+\lambda_{\max}(\tilde{U}^{T}+\tilde{U})$$

$$=H(\dots,\tilde{U}^{*})+n\lambda_{\max}(\tilde{U}^{T}+\tilde{U})+\operatorname{tr}(p_{\alpha}^{T}\tilde{U}\alpha)$$

$$=H(\dots,\tilde{U}^{*})+n\lambda_{\max}(\tilde{U}^{T}+\tilde{U})-2\operatorname{tr}(\tilde{U})$$

$$=H(\dots,\tilde{U}^{*})+n\lambda_{\max}(\tilde{U}^{T}+\tilde{U})-\operatorname{tr}(\tilde{U}^{T}+\tilde{U})$$

$$=H(\dots,\tilde{U}^{*})+\sum_{i=1}^{n}(\underbrace{\lambda_{\max}(\tilde{U}^{T}+\tilde{U})-\lambda_{i}(\tilde{U}^{T}+\tilde{U})}_{\geq 0}).$$

Thus, $H(\dots,\tilde{U}^{*})$ satisfies PMP.