Distributionally Robust Probabilistic Prediction
for Stochastic Dynamical Systems

A probabilistic prediction is a forecast that assigns probabilities to different possible outcomes, enabling better uncertainty quantification, risk assessment, and decision-making. Due to the ever-increasing demand for both accurate and reliable predictions against uncertainty [1], it has wide applications in fields such as epidemiology [2], climatology [3], and robotics [4], etc. The performance of a probabilistic prediction is evaluated by both sharpness and calibration of the predictive probability density function (pdf). The sharpness means how concentrated the predictive pdf is, and calibration concerns the statistical compatibility between the predictive pdf and the realized outcomes [5]. To simultaneously assess both sharpness and calibration, one usually assigns a numerical score to the predictive pdf and the realized value of the prediction target. If one has a further requirement that the expected scoring rule is maximized only when the predictive pdf equals the target’s real pdf almost everywhere, a strictly proper scoring rule (see in [6]) is needed.

The probabilistic prediction of a stochastic dynamical system (SDS) is to predict the conditional pdfs of future states based on some prior information about the system. One’s prior information is usually in the form of a nominal model, which is a simplified or idealized representation of the underlying system. While probabilistic prediction for a linear system with Gaussian noises is straightforward, it becomes computationally expensive once the nominal SDS is nonlinear or the noises are non-Gaussian. Therefore, a lot of research contributes to approximating the future pdf to balance the requirement of prediction precision and computation complexity [7].

Despite the significant achievements in the current prediction algorithms for SDSs, they are not guaranteed to perform well if the distributional information provided by a nominal model is inaccurate. It is ubiquitous for a nominal SDS to possess both inaccurate state evolution functions and incomplete distributional information of noises. For example, the nominal state evolution function may be a local linearization from a nonlinear model [8]. Additionally, one’s nominal information of the distributions of system noises is usually partial. Most of the time, only estimated mean and covariance are available in the nominal model, which is far from uniquely determining a pdf [9]. A nominal predictor can lead to misleading predictive pdfs even though the uncertainties on SDS are not prominent (as demonstrated by experiments in Sec. VIII). Addressing this issue is crucial for the widespread application of probabilistic prediction in model-uncertain scenarios. In this paper, the aim is to develop a framework to design probabilistic predictors for SDSs with worst-case performance guarantees.

To design probabilistic predictors against distributional uncertainties of model, we inherit the key notions from the distributionally robust optimization (DRO) [10]. DRO is a minimax-based mathematical framework for decision-making under uncertainty. It aims to find the decision that performs well across a range of possible pdfs in an ambiguity set. Nevertheless, what we need is to find the optimal probabilistic predictor that performs well over a set of SDSs. In DRPP, we generalize the idea of ambiguity sets to jointly describe one’s uncertainties of noises’ pdfs and state evolution functions.

Motivated by the aforementioned considerations, we propose a novel functional-maximin-based distributionally robust probabilistic prediction (DRPP) framework. In DRPP, the objective is the expected score under the worst-case SDS within an ambiguity set, and the predictor optimizes the objective over the function space of all pdfs on the state space. Compared to DRO, the proposed DRPP is significantly different in three perspectives. i) (Functional optimization) One needs to find the optimal pdfs over function spaces rather than optimizing decisions over Euclidean spaces. ii) (Conditional prediction) The predictor for an SDS is a map from the state space to the pdf space, rather than a single decision variable without conditioning on some context. iii) (Trajectory expectation) The objective is an expectation concerning all possible trajectories and control inputs, rather than an expectation over a single random variable. The main contributions of this paper are summarized as follows:

• •

  To the best of our knowledge, this is the first work on probabilistic prediction of SDSs that optimizes the worst-case performance over an ambiguity set. We have generalized the classical conic moment-based ambiguity set such that both the uncertainties of state evolution functions and system noises are quantified.

• •

  We develop a methodology to greatly reduce the complexity of solving DRPP. First, we use the principle of dynamic programming to prove that solving DRPP only needs to concern the one-step DRPP. Second, we exploit a necessary optimality condition to equivalently transform the maximin problem in function spaces into a convex-nonconvex minimax problem in Euclidean spaces.

• •

  We design two suboptimal DRPP predictors by relaxing constraints. Noise-DRPP is the optimal predictor when the nominal state evolution function is accurate, and Eig-DRPP is the optimal predictor when the predictor is constrained by an eigenvector restriction. Moreover, an explicit upper bound of the optimality gap between Eig-DRPP and the global optimal predictor is obtained.

• •

  We provide upper and lower bounds for the optimal worst-case value function of the DRPP. Utilizing these theoretical bounds, we synthesize practical rules to judge whether a pre-defined ambiguity set for SDS is too large or too small. Based on these rules, one can adjust the ambiguity set in real-time for better performance.

The rest of this article is organized as follows. Section II provides a brief review of the related work and Section III introduces preliminaries on probabilistic prediction and DRO. In Section IV, we formulate the proposed DRPP problem. In Section V, we introduce the main methodology to address the main challenges in the DRPP. Then, we relax the constraints to solve two suboptimal DRPPs with upper and lower bounds in Section VI and Section VII. Next, we conduct elaborate experiments and provide an application in Section VIII, followed by a conclusion in Section IX.

Suppose a random variable $\boldsymbol{x}$ takes value on $\mathcal{X}$ with distribution $p_{\boldsymbol{x}}$, the output of a probabilistic predictor is a predictive distribution $\hat{p}_{\boldsymbol{x}}\in\mathcal{P}(\mathcal{X})$, where $\mathcal{P}(\mathcal{X})$ denotes the set of all probability distributions over $\mathcal{X}$. After the value of $\boldsymbol{x}$ is materialized as $x$, a scoring rule,

assigns a numerical score $S(\hat{p}_{\boldsymbol{x}},x)$ to measure the quality of the predictive distribution $\hat{p}_{\boldsymbol{x}}$ on the realized value $x$. The expected scoring rule of $\mathcal{S}$, usually sharing the same operator $\mathcal{S}:\mathcal{P}(\mathcal{X})\times\mathcal{P}(\mathcal{X})\to\mathbb{R}$ but different operands, is

A scoring rule $\mathcal{S}$ is proper with respect to the distribution space $\mathcal{P}(\mathcal{X})$ if $\mathcal{S}(\hat{p}_{\boldsymbol{x}},p_{\boldsymbol{x}})\leq\mathcal{S}(p_{\boldsymbol{x}},p_{\boldsymbol{x}})$ holds for all $\hat{p}_{\boldsymbol{x}},p_{\boldsymbol{x}}\in\mathcal{P}(\mathcal{X})$. It is strictly proper if the equality holds only when $\hat{p}_{\boldsymbol{x}}$ equals $p_{\boldsymbol{x}}$ almost everywhere. For example, the logarithm score,

is one of the most celebrated scoring rules for being essentially the only local proper scoring rule up to equivalence [36].

Traditional stochastic optimization approaches assume a known probability distribution $P$ and solve problems

where $x\in\mathcal{X}$ represents the decision variable, $\xi$ is sampled from the distribution $P$, $h(x,\xi)$ is the objective function. However, in many real-world problems, the true probability distribution of the uncertain parameters is unknown or partially known. DRO addresses this issue by considering a family of distributions, known as an ambiguity set, rather than a single distribution. The goal is to find a solution that performs well for the worst-case distribution within this set. A general DRO problem is formulated as:

where $\mathcal{D}$ is the ambiguity set, a set of probability distributions that are close to the nominal distribution $P$ in some sense.

The ambiguity set $\mathcal{D}$ can be defined in various ways. One of the most popular ones is the moment-based sets, where distributions are constrained by their moments (e.g., mean, variance). For example, the celebrated conic moment-based ambiguity set is used in the seminal work [10] such that

where $\mathcal{X}$ is the nominal support, $\mu_{0}$ and $\Sigma_{0}$ are the nominal mean and covariance, $\gamma_{1},\gamma_{2}\in\mathbb{R}_{+}$ quantifies the uncertainties on the mean and covariance. Other popular ambiguity sets include structure-based sets, distance-based sets, etc. A thorough review of popular ambiguity sets is provided in [37, Sec. 3].

System model $(\mathcal{X},\mathcal{U},\Phi,\pi,T)$. Consider a discrete-time time-variant SDS with horizon $T\in\mathbb{Z}_{+}$, denoted as $\Phi$. The dynamics of $\Phi$ at time step $k\in\{0,\ldots,T-1\}$ is

where $\boldsymbol{x}_{k}$ is the system state taking values in $\mathcal{X}=\mathbb{R}^{d_{x}}$, the initial state is given as $x_{0}$, $f_{k}:\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{u}}\to\mathbb{R}^{d_{x}}$ is the state evolution function. $\boldsymbol{u}_{k}$ is the control input taking values in $\mathcal{U}=\mathbb{R}^{d_{u}}$, which is generated from a state-feedback control policy $\pi_{k}:\mathbb{R}^{d_{x}}\to\mathcal{P}(\mathbb{R}^{d_{u}})$, i.e., $\boldsymbol{u}_{k}\sim\pi_{k}(\cdot\mid x_{k})$. $\boldsymbol{w}_{k}$ is the system noise taking values in $\mathbb{R}^{d_{x}}$. Since the pair $(f_{k},p_{\boldsymbol{w}_{k}})$ determines $\Phi_{k}$, it is convenient to denote the pair $(f_{k},p_{\boldsymbol{w}_{k}})$ as $\Phi_{k}$. Without raising confusion, we also denote the state-control pair $(\boldsymbol{x}_{k},\boldsymbol{u}_{k})$ as $\boldsymbol{z}_{k}$.

Prediction model $(\mathcal{S},\mathscr{F})$. In this paper, we use the popular logarithm score to measure the performance of probabilistic predictions, i.e., $\mathcal{S}=\mathcal{L}$. Starting from the system’s initial state $x_{0}$, a probabilistic predictor keeps observing the states and the control inputs, aiming to predict the conditional pdfs of the next state. The probabilistic prediction of the SDS state at time step $k$ is a map $\mathscr{F}_{k}$ from the previous state and control input to a predictive conditional pdf for the next state, i.e.,

The prediction performance of $\mathscr{F}_{k}$ at time step $k$ is measured as $\mathcal{L}\left(\mathscr{F}_{k}(\boldsymbol{z}_{k}),\boldsymbol{x}_{k+1}\right)$. Next, we measure the prediction performance of a probabilistic predictor along an SDS state trajectory by the sum of one-step scores.

Because the logarithm scoring rule is strictly proper, the expected logarithm score can only attain its maximum when $\hat{p}_{\boldsymbol{w}}$ identifies with $p_{\boldsymbol{{w}}}$ almost everywhere. Therefore, a predictor needs full information on both the state evolution function and the noises’ distributions to maximize the expected scores, which is often impossible in practice. We detail this lack of information in the following assumption on the predictor’s prior information during the trajectory prediction.

The predictor’s prior information usually does not identify with the ground truth. To further quantify the uncertainty between the nominal model and the real model, we need to define an ambiguity set of SDSs.

The crucial spirit of an ambiguity set is, that although a nominal model rarely identifies with the ground truth, one can still have additional information or belief that the real model locates around the nominal model. For example, an upper bound for the operator norm $\|f_{k}-\bar{f}_{k}\|_{2}$ may be available when the nominal function is derived from a system identification procedure. Additionally, estimations and confidence regions (CRs) of the mean and covariance of system noise may be statistically inferred from historical data. Even if no supportive evidence exists to construct a reasonable ambiguity set, one can still subjectively select an ambiguity set.

Jointly considering the uncertainties of both state evolution functions and noises, we define the conic moment-based ambiguity set for SDSs.

When there is no confusion about the parameters that determine an ambiguity set, we often use $\mathcal{I}_{k}$ for simplification.

To begin with, we need to formulate the expected cumulative scores that a predictor can achieve from a given state-control pair. In reinforcement learning, the value function is a mathematical tool to estimate the expected cumulative reward for an agent following certain policy [38]. This concept can also be applied to estimate the performance of a probabilistic predictor. Starting from a time step $k\in\{0,\ldots,T-1\}$ with a state-control pair $z_{k}$, the expected performance of a predictor $\mathscr{F}$ is denoted by the value function $V^{\mathscr{F}}_{k}:\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{u}}\to\mathbb{R}$,

where the expectation is over the state-control trajectory. For convenience, we also use the notation

for the expected performance of $\mathscr{F}$ along the whole time horizon starting from initial state-control pair $z_{0}$.

Next, given an initial state-control pair $z_{0}$ and a sequence of ambiguity sets $\mathcal{I}_{0:T-1}$ of the underlying SDS $\Phi$, we are interested in designing probabilistic predictors that can perform well in the worst-case scenarios. To this end, a DRPP problem for SDS is formulated as a multistep functional maximin optimization:

where the inner minimization yields the worst-case performance for a given predictor over the ambiguity sets, the outer maximization essentially finds out the optimal predictive pdfs to maximize the worst-case performance. Please see Fig. 1 for the overall structure of this paper.

Given a state-control pair $z_{k}$ at time step $k\in\{0,\ldots,T-1\}$, we denote the optimal worst-case value function as

Leveraging the principle of dynamic programming, one can derive the following recursive equations for $V^{*}_{k}$,

Thanks to the principle of dynamic programming, one has that the optimal predictor at time step $k$ is

Conditioned on the previous state-control pair $z_{k}$, optimizing over the probabilistic predictor $\mathscr{F}_{k}$ is equivalent to optimizing over the predictive pdf. Therefore, conditioned on $z_{k}$ at time step $k\in\{0,\ldots,T-1\}$, we call the maximin optimization over the predictive pdf space $\mathcal{P}(\mathbb{R}^{d_{x}})$ and the ambiguity set $\mathcal{I}_{k}$ as a one-step DRPP:

where $\hat{p}_{k}(\cdot):=\mathscr{F}_{k}(z_{k})=\hat{p}_{\boldsymbol{x}_{k+1}\mid\boldsymbol{z}_{k}}(\cdot\mid z_{k})$ is an abbreviation for the predictive pdf of $\boldsymbol{x}_{k+1}$ conditioned on $z_{k}$. In the rest of this paper, we focus on solving the one-step DRPP $\mathbf{P}_{1}$.

Towards solving the one-step DRPP $\mathbf{P}_{1}$, the very first step is to canonicalize the conic moment-based ambiguity set $\mathcal{I}_{k}$ such that all the optimization variables appeared in the objective and constraints are consistent.

First, to canonicalize the first ambiguity constraint on $f_{k}$, given a state-control pair $z_{k}$ at time step $k$, we define the real one-step state evolution as $\nu_{k}=f_{k}(z_{k})$, and the nominal one-step state evolution as $\bar{\nu}_{k}=\bar{f}_{k}(z_{k})$. Then, this ambiguity constraint of $f_{k}$ is equivalent to describing the distance between $\nu_{k}$ and $\bar{\nu}_{k}$, i.e.,

Second, to canonicalize the remaining ambiguity constraints on $p_{\boldsymbol{w}_{k}\mid\boldsymbol{z}_{k}(\cdot\mid z_{k})}$, we need to keep them consistent with the pdf $p_{\boldsymbol{x}_{k+1}\mid\boldsymbol{z}_{k}}(\cdot\mid z_{k})$ in the optimization objective. For ease of notation, since the one-step DRPP problem is always conditioned on a given state-control pair $z_{k}$, we abbreviate the notations by denoting $p_{\boldsymbol{\xi}_{k}}(\cdot)=p_{\boldsymbol{x}_{k+1}\mid\boldsymbol{z}_{k}}(\cdot\mid z_{k})$ and $p_{\boldsymbol{w}_{k}}(\cdot)=p_{\boldsymbol{w}_{k}\mid\boldsymbol{z}_{k}}(\cdot\mid z_{k})$ without raising confusion. It follows that $p_{\boldsymbol{\xi}_{k}}(\xi)=p_{\boldsymbol{w}_{k}}(\xi-\nu_{k})$, and we shall align one with another. Here we align $p_{\boldsymbol{w}_{k}}$ with $p_{\boldsymbol{\xi}_{k}}$ and then optimize over $p_{\boldsymbol{\xi}_{k}}$. Finally, the inner minimization of the one-step DRPP is canonicalized as an infinite dimensional conic linear program:

Directly solving the worst-case pdf in the function space $\mathcal{P}(\mathbb{R}^{d_{x}})$ is difficult. As a first step towards dealing with this challenge, we leverage the dual analysis for problem $\mathbf{P}_{1}$. If a strong duality holds, one can optimize the Lagrange multipliers in the dual form. For the inner optimization (8), a classical conclusion is that $\bar{\Sigma}_{k}\succ 0$ is a sufficient condition for strong duality to hold [10, p. 5]. Let $r\in\mathbb{R},q\in\mathbb{R}^{d_{x}},Q_{i}\in S_{+}^{d_{x}},\kappa_{i}=\{P_{i}\in S_{+}^{d_{x}},p_{i}\in\mathbb{R}^{d_{x}},s_{i}\in\mathbb{R}\}$ for $i=1,2$ be the Lagrange multipliers, the dual form of (8) is

where $G(r,q,Q_{1},Q_{2},\kappa_{1},\kappa_{2})\!=\!-r\!+\!\bar{\mu}_{k}^{\top}\!(Q_{1}\!-\!Q_{2})\bar{\mu}_{k}-(\gamma_{2,k}Q_{1}-\gamma_{3,k}Q_{2}+P_{1})\cdot\bar{\Sigma}_{k}-P_{2}\cdot I+2p_{1}^{\top}\bar{\mu}_{k}-s_{1}\gamma_{1,k}+2p_{2}^{\top}\bar{\nu}_{k}-s_{2}\gamma_{0,k}\|z_{k}\|^{2}_{2}-(q+2p_{2})^{\top}\nu_{k}-\nu_{k}^{\top}(Q_{1}-Q_{2})\nu_{k}$. Please see Appendix B for a detailed derivation.

Unfortunately, we have to deal with an infinite number of constraints indexed by $w\in\mathbb{R}^{d_{x}}$ in problem $\mathbf{D}_{1}$. Although the dual problem cannot lead to a direct solution, it helps to develop an optimality condition for the probabilistic predictors in the following theorem.

Theorem 1 reduces the optimization complexity from the original function space $\mathcal{P}(\mathbb{R}^{d_{x}})$ to the Gaussian distribution family, more specifically, the parameter spaces of mean and covariance. Therefore, $\mathbf{D}_{1}$ can be transformed into a finite-dimensional optimization problem as follows

Please see Appendix D for a detailed derivation. However, this problem is now a nonconvex semidefinite optimization with respect to $\hat{\Sigma}_{k}$, solving a globally optimal solution is NP-hard in general. An alternative is to use numerical iterative algorithms to seek a locally optimal solution of (10). Nevertheless, the local optimality of the dual problem cannot guarantee the worst-case performance of $\mathbf{P}_{1}$. To deal with this hardness, we avoid directly solving the dual problem by applying the optimality condition in Theorem 1 to $\mathbf{P}_{1}$.

Substituting the parameterized $\hat{p}_{k}$ into the primal one-step DRPP $\mathbf{P}_{1}$, we have the objective being as

Next, notice that $\mathbb{E}_{\xi\sim p_{\boldsymbol{\xi}_{k}}}[\xi]=\mu_{k}+\nu_{k}$ and $\mathbb{E}_{\xi\sim p_{\boldsymbol{\xi}_{k}}}(\xi-\nu_{k}\!-\mu_{k})(\xi-\nu_{k}-\mu_{k})^{\top}=\Sigma_{k}$, one has $\mathbb{E}_{\xi\sim p_{\boldsymbol{\xi}_{k}}}(\xi-\hat{\mu}_{k})\hat{\Sigma}_{k}^{-1}(\xi-\hat{\mu}_{k})^{\top}$ equals to $[\Sigma_{k}+(\mu_{k}+\nu_{k}-\hat{\mu}_{k})(\mu_{k}+\nu_{k}-\hat{\mu}_{k})^{\top}]\cdot\hat{\Sigma}_{k}^{-1}$. The functional minimax $\mathbf{P}_{1}$ is then equivalent to a minimax optimization in Euclidean spaces as follows

Substituting the expressions of $\Sigma_{k}^{*}$ and $\hat{\mu}_{k}^{*}$ into $\mathbf{P}_{2}$, it is equivalent to a convex-nonconcave minimax optimization:

where $a=\mu_{k}-\bar{\mu}_{k},b=\nu_{k}-\bar{\nu}_{k}$. Notice that the inner maximization is an indefinite quadratic constrained quadratic programming (QCQP) problem, which is still short of efficient algorithms to solve this kind of minimax optimization.

Given any $\hat{\Sigma}_{k}\succ 0$ and $b\neq 0$, one can first solve $a$ by concerning the following quadratic constrained linear programming (QCLP) problem:

Substituting the optimal $a^{*}$ into the problem, the inner optimization of (11) can be further reformulated as a convex maximization problem:

where $A=\hat{\Sigma}_{k}^{-1}\bar{\Sigma}_{k}\hat{\Sigma}_{k}^{-1}$, $B=\hat{\Sigma}_{k}^{-1}$, and $\alpha=2\sqrt{\gamma_{1,k}}$.

As we have shown, solving the one-step DRPP $\mathbf{P}_{1}$ is equivalent to finding a global minimax point (Stackelberg equilibrium) for a convex-nonconcave minimax problem where the outer minimization is a semidefinite programming and the inner problem is to maximize a non-smooth convex function subject to quadratic constraint. It is known that convex maximization is NP-hard in very simple cases such as quadratic maximization over a hypercube, and even verifying local optimality is NP-hard [40]. Thus finding a global minimax point of a one-step DRPP problem is also NP-hard.

Even worse, finding a local minimax point of a general constrained minimax optimization is of fundamental hardness. Not only the local minimax point may not exist [41], but any gradient-based algorithm needs exponentially many queries in the dimension and $\epsilon^{-1}$ to compute an $\epsilon$-approximate first-order local stationary point [42]. Finally, the local stationary point is not always guaranteed to be a local minimax point [43, 44].

In summary, globally finding a minimax point of a one-step DRPP is computationally intractable, and gradient-based algorithms aiming for local minimax points are limited by both convergence speed and performance guarantee. Therefore, a reasonable goal for the one-step DRPP is to derive upper or lower bounds of the optimal worst-case value function $V_{k}^{*}$ by relaxing the constraints imposed on the predictors or SDSs.

If we shrink the ambiguity set of one-step DRPP $\mathbf{P}_{1}$ by assuming that the nominal one-step state evolution $\bar{\nu}_{k}$ is equal with the nominal one-step state evolution $\nu_{k}$, there is $w_{k}=x_{k+1}-\bar{\nu}_{k},$ which means that the value of $w_{k}$ is precisely known after $x_{k+1}$ been observed.

Conditioned on $z_{k}$, we let $\mu_{k}$ be the conditional expectation of $\boldsymbol{w}_{k}$ and abbreviate the conditional pdf $p_{\boldsymbol{w}_{k}\mid\boldsymbol{z}_{k}}(\cdot\mid z_{k})$ as $p_{\boldsymbol{w}_{k}}$. Different from the general condition as discussed in Remark 3 where one should not formulate the inner minimization over $p_{\boldsymbol{w}_{k}}$, now the inner minimization can be structured as a canonical moment problem on $\mu_{k}$ and $p_{\boldsymbol{w}_{k}}$:

With the outer maximization imposed on $\hat{p}_{k}$, we can get a reformulated one-step DRPP whose optimal worst-case value function is an upper bound of $V^{*}(z_{k})$.

Taking the dual form of the inner moment problem and joining it with the outer maximization, we transform the primal maximin problem into the following maximization problem:

where $G(r,q,Q_{1},Q_{2},P,p,s)=-r+\bar{\mu}_{k}^{\top}(Q_{1}-Q_{2})\bar{\mu}_{k}-(\gamma_{2,k}Q_{1}-\gamma_{3,k}Q_{2}+P)\cdot\bar{\Sigma}_{k}+2p^{\top}\bar{\mu}_{k}-s\gamma_{1,k}$. Please see Appendix F for a detailed derivation for the dual problem.