Distributionally Robust Optimization

Abstract

This chapter presents a class of distributionally robust optimization problems in which a decision-maker has to choose an action in an uncertain environment. The decision-maker has a continuous action space and aims to learn her optimal strategy. The true distribution of the uncertainty is unknown to the decision-maker. This chapter provides alternative ways to select a distribution based on empirical observations of the decision-maker. This leads to a distributionally robust optimization problem. Simple algorithms, whose dynamics are inspired from the gradient flows, are proposed to find local optima. The method is extended to a class of optimization problems with orthogonal constraints and coupled constraints over the simplex set and polytopes. The designed dynamics do not use the projection operator and are able to satisfy both upper- and lower-bound constraints. The convergence rate of the algorithm to generalized evolutionarily stable strategy is derived using a mean regret estimate. Illustrative examples are provided.

1. Introduction

Robust optimization can be defined as the process of determining the best or most effective result, utilizing a quantitative measurement system under worst case uncertain functions or parameters. The optimization may occur in terms of best robust design, net cash flows, profits, costs, benefit/cost ratio, quality-of-experience, satisfaction, end-to-end delay, completion time, etc. Other measurement units may be used, such as units of production or production time, and optimization may occur in terms of maximizing production units, minimizing processing time, production time, maximizing profits, or minimizing costs under uncertain parameters. There are numerous techniques of robust optimization methods such as robust linear programming, robust dynamic programming, robust geometric programming, queuing theory, risk analysis, etc. One of the main drawbacks of the robust optimization is that the worst scenario may be too conservative. The bounds provided by the worst case scenarios may not be useful in many interesting problems (see the wireless communication example provided below). However, distributionally robust optimization is not based on the worst case parameters. The distributional robustness method is based the probability distribution instead of worst parameters. The worse case distribution within a certain carefully designed distributional uncertainty set may provide interesting features. Distributionally robust programming can be used not only to provide a distributionally robust solution to a problem when the true distribution is unknown, but it also can, in many instances, give a general solution taking into account some risk. The presented methodology is simple and reduces significantly the dimensionality of the distributionally robust optimization. We hope that the designs of distributionally robust programming presented here can help designers, engineers, cost–benefit analyst, managers to solve concrete problems under unknown distribution.

The rest of the chapter is organized as follows. Section 2 presents some preliminary concepts of distributionally robust optimization. A class of constrained distributionally robust optimization problems are presented in Section 3. Section 4 focuses on distributed distributionally robust optimization. Afterwards, illustrative examples in distributed power networks and in wireless communication networks are provided to evaluate the performance of the method. Finally, prior works and concluding remarks are drawn in Section 5.

Notation: Let $\mathbb{R}$, ${\mathbb{R}}_{+}$, denote the set of real and non-negative real numbers, respectively, $\left(\mathrm{\Omega },d\right)$ be a separable completely metrizable topological space with $d:\mathrm{\Omega }\times \mathrm{\Omega }\to {\mathbb{R}}_{+}$ a metric (distance). Let $\mathcal{P}\left(\mathrm{\Omega }\right)$ be the set of all probability measures over $\mathrm{\Omega }.$

2. Distributionally robust optimization

This section introduces distributionally robust optimization models. We will first present a generic formulation of the problem. Then, individual components of the optimization and their solvability issues via equivalent formulations will be discussed.

2.1. Model

Consider a decision-maker who wants to select an action $a\in \mathcal{A}\subset {\mathbb{R}}^n$ in order to optimize her objective $r\left(a,\omega \right),$ where $\omega$ is an uncertain parameter. The information structure is the following:

• The true distribution of $\omega$ is not known to the decision-maker.

• The upper/lower bound (if any) of $\omega$ are unknown to the decision-maker.

• The decision-maker can measure/observe realization of the random variable $\omega .$

The decision-maker chooses to experiment several trials and obtains statistical realizations of $\omega$ from measurements. The measurement data can be noisy, imperfect and erroneous. Then, an empirical distribution (or histogram) $m$ is built from the realizations of $\omega .$ However, $m$ is not the true distribution of the random variable $\omega ,$ and $m$ may not be a reliable measure due to statistical, bias, measurement, observation or computational errors. Therefore, the decision-maker is facing a risk. The risk-sensitive decision-maker should decide action that improves the performance of ${\mathbb{E}}_{\tilde{m}}r\left(a,\omega \right)$ among alternative distributions $\tilde{m}$ within a certain level of deviation $\rho >0$ from the distribution $m.$ The distributionally robust optimization problem is therefore formulated as

$${sup}_{a\in \mathcal{A}}{inf}_{\tilde{m}\in B_{\rho }\left(m\right)}{\mathbb{E}}_{\omega \sim \tilde{m}}r\left(a,\omega \right).$$E1

where $B_{\rho }\left(m\right)$ is the uncertainty set of alternative admissible distributions from $m$ within a certain radius $\rho >0.$ Different distributional uncertainty sets are presented: the $f$-divergence and the Wasserstein metric, defined below.

2.1.1. $f$-divergence

We introduce the notion of $f-$ divergence which will be used to compute the discrepancy between probability distributions.

Definition 1. Let $m$ and $\tilde{m}$ be two probability measures over $\mathrm{\Omega }$ such that $m$ is absolutely continuous with respect to $\tilde{m}.$ Let $f$ be a convex function. Then, the $f$-divergence between $m$ and $\tilde{m}$ is defined as follows:

$$D_f\left(m\parallel \tilde{m}\right)\equiv {\int }_{\mathrm{\Omega }}f\left(\frac{\mathit{dm}}{d\tilde{m}}\right)d\tilde{m}-f\left(1\right),$$

where $\frac{\mathit{dm}}{d\tilde{m}}$ is the Radon-Nikodym derivative of the measure $m$ with the respect the measure $\tilde{m}.$

By Jensen’s inequality:

$$\begin{array}{lll}{D}_f\left(m\parallel \tilde{m}\right)& =& {\int }_{\mathrm{\Omega }}f\left(\frac{\mathit{dm}}{d\tilde{m}}\right)d\tilde{m}-f\left(1\right)\\ & \ge & f\left({\int }_{\mathrm{\Omega }}\frac{\mathit{dm}}{d\tilde{m}}d\tilde{m}\right)-f\left(1\right)\\ & =& f\left({\int }_{\mathrm{\Omega }}\mathit{dm}\right)-f\left(1\right)\\ & =& f\left(1\right)-f\left(1\right)=0.\end{array}$$E2

Thus, $D_f\left(m\parallel \tilde{m}\right)\ge 0$ for any convex function $f.$ Note however that, the $f-$ divergence $D_f\left(m\parallel \tilde{m}\right)$ is not a distance (for example, it does not satisfy the symmetry property). Here the distributional uncertainty set imposed to the alternative distribution $\tilde{m}$ is given by

$$B_{\rho }\left(m\right)=\left\{\tilde{m}|\tilde{m}\left(.\right)\ge 0,{\int }_{\mathrm{\Omega }}d\tilde{m}=\tilde{m}\left(\mathrm{\Omega }\right)=1,D_f\left(\tilde{m}\Vert m\right)\le \rho \right\}.$$

Example 1. From the notion of $f-$ divergence one can derive the following important concept:

• $\alpha$-divergence for

$$f\left(a\right)=\left\{\begin{array}{ll}\frac{4}{\left(\alpha +1\right)\left(1-\alpha \right)}\left(1-a^{\frac{\alpha +1}{2}}\right)& \mathit{\text{if}}\alpha \notin \left\{-1,+1\right\},\\ aloga& \mathit{\text{if}}\alpha =1,\\ -loga& \mathit{\text{if}}\alpha =-1,\end{array}\right.$$

• In particular, Kullback–Leibler divergence (or relative entropy) is retrieved as $\alpha$ goes to $1.$

2.1.2. Wasserstein metric

The Wasserstein metric between two probability distributions $\tilde{m}$ and $m$ is defined as follows:

Definition 2. For $m,\tilde{m}\in \mathcal{P}\left(\mathrm{\Omega }\right),$ let $\mathrm{\Pi }\left(\tilde{m};m\right)$ be the set of all couplings between $m$ and $\tilde{m}.$ That is,

$$\left\{\pi \in \mathcal{P}\left(\mathrm{\Omega }\times \mathrm{\Omega }\right),|\pi \left(A\times \mathrm{\Omega }\right)=m\left(A\right),\pi \left(\mathrm{\Omega }\times B\right)=\tilde{m}\left(B\right),\left(A,B\right)\in {\mathcal{B}}^2\left(\mathrm{\Omega }\right)\right\}.$$

$\mathcal{B}\left(\mathrm{\Omega }\right)$ denotes the measurable sets of $\mathrm{\Omega }.$ Let $\theta \in \left[1,\infty \right].$ The Wasserstein metric between $\tilde{m}$ and $m$ is defined as

$$W_{\theta }\left(\tilde{m};m\right)=\underset{\pi \in \mathrm{\Pi }\left(\tilde{m};m\right)}{inf}\parallel d{\parallel }_{L_{\pi }^{\theta }}=\underset{\pi \in \mathrm{\Pi }\left(\tilde{m};m\right)}{inf}{\int }_{\left(a,b\right)}{d}^{\theta }\left(a,b\right)\pi \left(\mathit{da},\mathit{db}\right),$$

It is well-known that for every $\theta \ge 1,$ $W_{\theta }\left(\tilde{m};m\right)$ is a true distance in the sense that it satisfies the following three axioms:

• positive-definiteness,

• the symmetry property,

• the triangle inequality.

Note that $\tilde{m}$ is not necessarily absolutely continuous with respect to $m.$ Now the distributional uncertainty/constraint set is the set of all possible probability distributions within a $L^{\theta }$-Wasserstein distance below $\rho .$

$${\tilde{B}}_{\rho }\left(m\right)=\left\{\tilde{m}|{\int }_{\mathrm{\Omega }}d\tilde{m}=\tilde{m}\left(\mathrm{\Omega }\right)=1,W_{\theta }\left(\tilde{m};m\right)\le \rho \right\},$$

Note that, if $m$ is a random measure (obtained from a sampled realization), we use the expected value of the Wasserstein metric.

Example 2. The $L^{\theta }$-Wasserstein distance between two Dirac measures ${\delta }_{{\omega }_0}$ and ${\delta }_{{\tilde{\omega }}_0}$ is $W_{\theta }\left({\delta }_{{\omega }_0},{\delta }_{{\tilde{\omega }}_0}\right)=d\left({\omega }_0,{\tilde{\omega }}_o\right).$ More generally, for $K\ge 2,$ the $L^2$-Wasserstein distance between empirical measures ${\mu }_K=\frac{1}{K}{\sum }_{k=1}^K{\delta }_{{\omega }_k}$ and ${\nu }_K=\frac{1}{K}{\sum }_{k=1}^K{\delta }_{{\tilde{\omega }}_k}$ is $W_2^2\left({\mu }_K,{\nu }_K\right)\le \frac{1}{K}{\sum }_{i=1}^K{\left[{\omega }_k-{\tilde{\omega }}_k\right]}^2.$

We have defined $B_{\rho }\left(m\right)$ and ${\tilde{B}}_{\rho }\left(m\right).$ The goal now is to solve (1) under both $f-$ divergence and Wasserstein metric. One of the difficulties of the problem is the curse of dimensionality. The distributionally robust optimization problem (1) of the decision-maker is an infinite-dimensional robust optimization problem because $B_{\rho }$ is of infinite dimensions. Below we will show that (1) can be transformed into an optimization in the form of $supinfsup.$ The latter problem has three alternating terms. Solving this problem requires a triality theory.

2.2. Triality theory

We first present the duality gap and develop a triality theory to solve equivalent formulations of (1). Consider uncoupled domains ${\mathcal{A}}_i,i\in \left\{\mathrm{1,2,3}\right\}.$ For a general function $r_2,$ one has

$$\underset{a_2\in {\mathcal{A}}_2}{sup}\underset{a_1\in {\mathcal{A}}_1}{inf}{r}_2\left(a_1,a_2\right)\le \underset{a_1\in {\mathcal{A}}_1}{inf}\underset{a_2\in {\mathcal{A}}_2}{sup}{r}_2\left(a_1,a_2\right)$$

and the difference

$$\underset{a_1\in {\mathcal{A}}_1}{min}\underset{a_2\in {\mathcal{A}}_2}{max}{r}_2\left(a_1,a_2\right)-\underset{a_2\in {\mathcal{A}}_2}{max}\underset{a_1\in {\mathcal{A}}_1}{min}{r}_2\left(a_1,a_2\right),$$

is called duality gap. As it is widely known in duality theory from Sion’s Theorem [1] (which is an extension of von Neumann minimax Theorem) the duality gap vanishes, for example for convex-concave function, and the value is achieved by a saddle point in the case of non-empty convex compact domain.

Triality theory focuses on optimization problems of the forms: $supinfsup$ or $infsupinf.$ The term triality is used here because there are three key alternating terms in these optimizations.

Proposition 1. Let $\left(a_1,a_2,a_3\right)\mapsto r_3\left(a_1,a_2,a_3\right)\in \mathbb{R}$ be a function defined on the product space ${\prod }_{i=1}^3{\mathcal{A}}_i.$ Then, the following inequalities hold:

$$\begin{array}{l}{sup}_{a_2\in {\mathcal{A}}_2}{inf}_{a_1\in {\mathcal{A}}_1,a_3\in {\mathcal{A}}_3}{r}_3\left(a_1,a_2,a_3\right)\le \\ {inf}_{a_3\in {\mathcal{A}}_3}{sup}_{a_2\in {\mathcal{A}}_2}{inf}_{a_1\in {\mathcal{A}}_1}{r}_3\left(a_1,a_2,a_3\right)\le \\ {inf}_{a_1\in {\mathcal{A}}_1,a_3\in {\mathcal{A}}_3}{sup}_{a_2\in {\mathcal{A}}_2}{r}_3\left(a_1,a_2,a_3\right),\end{array}$$E3

and similarly

$$\begin{array}{l}{sup}_{a_1\in {\mathcal{A}}_1,a_3\in {\mathcal{A}}_3}{inf}_{a_2\in {\mathcal{A}}_2}{r}_3\left(a_1,a_2,a_3\right)\le \\ {sup}_{a_3\in {\mathcal{A}}_3}{inf}_{a_2\in {\mathcal{A}}_2}{sup}_{a_1\in {\mathcal{A}}_1}{r}_3\left(a_1,a_2,a_3\right)\le \\ {inf}_{a_2\in {\mathcal{A}}_2}{sup}_{a_1\in {\mathcal{A}}_1,a_3\in {\mathcal{A}}_3}{r}_3\left(a_1,a_2,a_3\right).\end{array}$$E4

Proof. Define

$$\widehat{g}\left(a_2,a_3\right)≔\underset{a_1\in {\mathcal{A}}_1}{inf}{r}_3\left(a_1,a_2,a_3\right).$$

Thus, for all $a_2,a_3,$ one has $\widehat{g}\left(a_2,a_3\right)\le r_3\left(a_1,a_2,a_3\right).$ It follows that, for any $a_1,a_3,$

$$\underset{a_2\in {\mathcal{A}}_2}{sup}\widehat{g}\left(a_2,a_3\right)\le \underset{a_2\in {\mathcal{A}}_2}{sup}{r}_3\left(a_1,a_2,a_3\right).$$

Using the definition of $\widehat{g},$ one obtains

$$\underset{a_2\in {\mathcal{A}}_2}{sup}\underset{a_1\in {\mathcal{A}}_1}{inf}{r}_3\left(a_1,a_2,a_3\right)\le \underset{a_2\in {\mathcal{A}}_2}{sup}{r}_3\left(a_1,a_2,a_3\right),\forall a_1,a_3.$$

Taking the infimum in $a_1$ yields:

$$\underset{a_2\in {\mathcal{A}}_2}{sup}\underset{a_1\in {\mathcal{A}}_1}{inf}{r}_3\left(a_1,a_2,a_3\right)\le \underset{a_1\in {\mathcal{A}}_1}{inf}\underset{a_2\in {\mathcal{A}}_2}{sup}{r}_3\left(a_1,a_2,a_3\right),\forall a_3$$E5

Now, we use two operations for the variable $a_3$:

• Taking the infimum in the inequality (5) in $a_3$ yields

$$\underset{a_3\in {\mathcal{A}}_3}{inf}\underset{a_2\in {\mathcal{A}}_2}{sup}\underset{a_1\in {\mathcal{A}}_1}{inf}{r}_3\left(a_1,a_2,a_3\right)\le \underset{a_3\in {\mathcal{A}}_3}{inf}\underset{a_1\in {\mathcal{A}}_1}{inf}\underset{a_2\in {\mathcal{A}}_2}{sup}{r}_3\left(a_1,a_2,a_3\right)$$

$$=\underset{\left(a_1,a_3\right)\in {\mathcal{A}}_1\times {\mathcal{A}}_3}{inf}\underset{a_2\in {\mathcal{A}}_2}{sup}{r}_3\left(a_1,a_2,a_3\right),$$

which proves the second part of the inequalities (3). The first part of the inequalities (3) follows immediately from (5).

• Taking the supremum in inequality (5) in $a_3$ yields

$$\underset{\left(a_2,a_3\right)\in {\mathcal{A}}_2\times {\mathcal{A}}_3}{sup}\underset{a_1\in {\mathcal{A}}_1}{inf}{r}_3\left(a_1,a_2,a_3\right)\le \underset{a_3\in {\mathcal{A}}_3}{sup}\underset{a_1\in {\mathcal{A}}_1}{inf}\underset{a_2\in {\mathcal{A}}_2}{sup}{r}_3\left(a_1,a_2,a_3\right),$$

which proves the first part of the inequalities (4). The second part of the inequalities (4) follows immediately from (5).

This completes the proof.

2.3. Equivalent formulations

Below we explain how the dimensionality of problem (1) can be significantly reduced using a representation by means of the triality theory inequalities of Proposition 1.

2.3.1. $f$-divergence

Interestingly, the distributionally robust optimization problem (1) under $f$-divergence is equivalent to the finite dimensional stochastic optimization problem (when $\mathcal{A}$ are of finite dimensions). To see this, the original problem need to be transformed. Let us introduce the likelihood functional $L\left(\tilde{\omega }\right)=\frac{d\tilde{m}}{\mathit{dm}}\left(\tilde{\omega }\right),$ and set

$$L_{\rho }\left(m\right)=\left\{L|{\int }_{\tilde{\omega }}f\left(L\left(\tilde{\omega }\right)\right)\mathit{dm}-f\left(1\right)\le \rho ,{\int }_{\tilde{\omega }}L\left(\tilde{\omega }\right)\mathit{dm}\left(\tilde{\omega }\right)=1\right\}.$$

Then, the Lagrangian of the problem is

$$\begin{array}{ll}\tilde{r}\left(a,L,\lambda ,\mu \right)=& {\int }_{\tilde{\omega }}r\left(a,\tilde{\omega }\right)L\left(\tilde{\omega }\right)\mathit{dm}\left(\tilde{\omega }\right)\\ & -\lambda \left(\rho +f\left(1\right)-{\int }_{\tilde{\omega }}f\left(L\left(\tilde{\omega }\right)\right)\mathit{dm}\left(\tilde{\omega }\right)\right)\\ & -\mu \left(1-{\int }_{\tilde{\omega }}L\left(\tilde{\omega }\right)\mathit{dm}\left(\tilde{\omega }\right)\right),\end{array}$$

where $\lambda \ge 0$ and $\mu \in \mathbb{R}.$ The problem becomes

$$\left\{{sup}_a\right.{inf}_{L\in L_{\rho }\left(m\right)}{sup}_{\lambda \ge 0,\mu \in \mathbb{R}}\tilde{r}\left(a,L,\lambda ,\mu \right).$$E6

A full understanding of problem $\left(6\right)$ requires a triality theory (not a duality theory). The use of triality theory leads to the following equation:

$$\left\{{sup}_{a\in \mathcal{A}}\right.{inf}_{\tilde{m}\in B_{\rho }\left(m\right)}{\mathbb{E}}_{\tilde{m}}\left[r\right]={sup}_{a\in \mathcal{A},\lambda \ge 0,\mu \in \mathbb{R}}{\mathbb{E}}_{m}h,$$E7

where $h$ is the integrand function $-\lambda \left(\rho +f\left(1\right)\right)-\mu -\lambda f^{\ast }\left(\frac{r+\mu }{-\lambda }\right),$ where $f^{\ast }$ is Legendre-Fenchel transform of $f$ defined by

$$f^{\ast }\left(\xi \right)=\underset{L}{sup}\left[⟨L,\xi ⟩-f\left(L\right)\right]=-\underset{L}{inf}\left[f\left(L\right)-⟨L,\xi ⟩\right].$$E8

Note that the righthand side of (7) is of dimension $n+2,$ which reduces considerably the dimensionality of the original problem (1).

2.3.2. Wasserstein metric

Similarly, the distributionally robust optimization problem under Wasserstein metric is equivalent to the finite dimensional stochastic optimization problem (when $\mathcal{A}$ is a set of finite dimension). If the function $\omega \mapsto r\left(a,\omega \right)$ is upper semi-continuous and $\left(\mathrm{\Omega },d\right)$ is a Polish space then the Wasserstein distributionally robust optimization problem is equivalent to

$$\left\{\begin{array}{l}{sup}_{a\in \mathcal{A}}{inf}_{\tilde{m}\in {\tilde{B}}_{\rho }\left(m\right)}{\mathbb{E}}_{\tilde{m}}\left[r\right]={sup}_{a\in \mathcal{A}}{sup}_{\lambda \ge 0}{\mathbb{E}}_m\left[\tilde{h}\right],\\ \tilde{h}={\mathit{\lambda \rho }}^{\theta }+\mu +{sup}_{\widehat{\omega }\in \mathrm{\Omega }}\left[r\left(a,\omega \right)-\mu -\lambda d^{\theta }\left(\omega ,\widehat{\omega }\right)\right];\end{array}\right.$$E9

The next subsection presents algorithms for computing a distributionally robust solution from the equivalent formulations above.

2.4. Learning algorithms

Learning algorithms are crucial for finding approximate solutions to optimization and control problems. They are widely used for seeking roots/kernel of a function and for finding feasible solutions to variational inequalities. Practically, a learning algorithm generates a certain trajectory (or a set of trajectories) toward a potential approximate solution. Selecting a learning algorithm that has specific properties such as better accuracy, more stability, less-oscillatory and quick convergence is a challenging task [2, 3, 4, 5]. From the calculus of variations point of view, however, a learning algorithm generates curves. Therefore, selecting an algorithm among the others leads to an optimal control problem on the spaces of curves. Hence, it is natural to use optimal control theory to derive faster algorithms for a family of curves. Bergman-based algorithms and risk-aware version of it are introduced below to meet specific properties. We start by introducing the Bregman divergence.

Definition 3. The Bregman divergence $d_g:\mathcal{A}\times \mathcal{A}\to \mathbb{R}$ is defined on a differentiable strictly convex function $g:\mathcal{A}\to \mathbb{R}.$ For two points $\left(a,b\right)\in {\mathcal{A}}^2$, it measures the gap between $g\left(a\right)$ and the first-order Taylor expansion of g around $a$ evaluated at $b$

$$d_g\left(a,b\right)≔g\left(a\right)-g\left(b\right)-⟨\nabla g\left(b\right),a-b⟩.$$

Example 3. From the Bregman divergence one gets other features by choosing specific functions $g:$

• If $g\left(a\right)={\sum }_{i=1}^n{a}_i^2$ then the Bregman divergence $d_g\left(a,b\right)={\sum }_{i=1}^n{\left(a_i-b_i\right)}^2$ is the squared standard Euclidean distance.

• If $g\left(a\right)={\sum }_{i=1}^n{a}_{i}log{a}_i$ is defined on the relative interior of the simplex, i.e., $a\in \left\{b|b\in {\left(0,1\right)}^n,{\sum }_{i=1}^n{b}_i=1\right\}$ then the Bregman divergence $d_g\left(a,b\right)={\sum }_{i=1}^n{a}_{i}log\left(\frac{a_i}{b_i}\right),$ is the Kullback–Leibler divergence.

We are now ready to define algorithms for solving the righthand side of (7) and (9). One of the key approaches for error quantification of the algorithm with respect to the distributionally robust optimum is the so-called average regret. When the regret vanishes one gets close to a distributionally robust optimum.

Definition 4. The average regret of an algorithm which generates the trajectory $a\left(t\right)=\left(\tilde{a}\left(t\right),\lambda \left(t\right),\mu \left(t\right)\right)$ within $\left[t_0,T\right],t_0>0$ is

$${\mathit{\text{regret}}}_T≔\frac{1}{T-t_0}{\int }_{t_0}^T\left[\underset{b\in \mathcal{A}\times {\mathbb{R}}_{+}\times \mathbb{R}}{max}{\mathbb{E}}_{m}h\left(b,\omega \right)\right]-{\mathbb{E}}_{m}h\left(a\left(t\right),\omega \right)\mathit{dt}$$

2.4.1. Armijo gradient flow

Algorithm 1. The Armijo’s gradient pseudocode is as follows:

1: Procedure Armijo gradient $\left(a\left(0\right),\mathrm{ϵ},T,g,m,h\right)⊳$ The Armijo’s gradient starting from $a\left(0\right)$ within $\left[0,T\right]$

2:       $a\leftarrow a\left(0\right)$

3:       while $\mathit{\text{regret}}>\mathrm{ϵ}$ and $t\le T$ do $⊳$ We have the answer if regret is 0

4:              Compute $a\left(t\right)$ solution of (10)

5:              Compute ${\mathit{\text{regret}}}_t$

6:       end while

7:       return $a\left(t\right),{\mathit{\text{regret}}}_t$ $⊳$ get a(t) and the regret

8: end procedure

Proposition 2. Let $a\mapsto {\mathbb{E}}_{m}h\left(a,\omega \right):{\mathbb{R}}^{n+2}\to \mathbb{R}$ be a concave function that has a unique global maximizer $a^{\ast }.$ Assume that $a^{\ast }$ be a feasible action profile, i.e., $a^{\ast }\in \mathcal{A}.$ Consider the continuous time analogue of the Armijo gradient flow [6], which is given by

$$\begin{array}{ll}\frac{d}{\mathit{dt}}a\left(t\right)={\left[{\nabla }^{2}g\right]}^{-1}.{\nabla }_a{\mathbb{E}}_{m}h\left(a\left(t\right),\omega \right),& \\ a\left(0\right)=a_0\in {\mathbb{R}}^{n+2},& \end{array}$$E10

where $a\left(0\right)=a_0$ is the initial point of the algorithm and $g$ is a strictly convex function on $a.$ Let $a\left(t\right)$ be the solution to (10).

Then the average regret within $\left[t_0,T\right],t_0>0$ is bounded above by

$${\mathit{\text{regret}}}_T≔\frac{1}{T-t_0}{\int }_{t_0}^T{\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a\left(t\right),\omega \right)\right]\mathit{dt}\le d_g\left(a^{\ast },a_0\right)\frac{log\frac{T}{t_0}}{T-t_0}.$$

Proof. Let

$$W\left(a\left(t\right)\right)=t{\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a\left(t\right),\omega \right)\right]+d_g\left(a^{\ast },a\left(t\right)\right),$$

where $a$ is solution to (10). The function $W$ is positive and $\frac{d}{\mathit{dt}}W={\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a\left(t\right),\omega \right)\right]-t⟨{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right),g_{\mathit{aa}}^{-1}{\mathbb{E}}_m{\nabla }_{a}h\left(a\left(t\right),\omega \right)⟩+\frac{d}{\mathit{dt}}{d}_g\left(a^{\ast },a\left(t\right)\right).$ By concavity of ${\mathbb{E}}_{m}h\left(a,\omega \right)$ one has

$$⟨{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right),\left(a^{\ast }-a\right)⟩\ge {\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a,\omega \right)\right],\forall a.$$

On the other hand,

$$\begin{array}{ll}\frac{d}{\mathit{dt}}{d}_g\left(a^{\ast },a\left(t\right)\right)=-\dot{a}{g}_a\left(a\right)-⟨g_{\mathit{aa}}\dot{a},a-a^{\ast }⟩+g_a\dot{a}& \\ =-⟨g_{\mathit{aa}}\dot{a},a-a^{\ast }⟩=-⟨{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right),a^{\ast }-a⟩.& \end{array}$$E11

Hence,

$$\begin{array}{ll}\frac{d}{\mathit{dt}}W\le ⟨{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right),\left(a^{\ast }-a\right)⟩& \\ -t⟨{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right),g_{\mathit{aa}}^{-1}{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right)⟩& \\ -⟨{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right),a^{\ast }-a⟩& \\ =-t⟨{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right),g_{\mathit{aa}}^{-1}{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right)⟩\le 0,& \end{array}$$E12

where the last inequality is by convexity of $g.$ It follows that $\frac{d}{\mathit{dt}}W\left(a\left(t\right)\right)\le 0$ along the path of the gradient flow. This decreasing property implies $0\le W\left(a\left(t\right)\right)\le W\left(a\left(0\right)\right)=d_g\left(a^{\ast },a\left(0\right)\right).$ In particular, $0\le t{\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a,\omega \right)\right]\le W\left(a\left(0\right)\right)<+\infty .$ Thus, the error to the value ${\mathbb{E}}_{m}h\left(a^{\ast },\omega \right)$ is bounded by

$$0\le {\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a,\omega \right)\right]\le \frac{W\left(a\left(0\right)\right)}{t}.$$

The announced result on the regret follows by integration over $\left[t_0,T\right]$ and by averaging. This completes the proof.

Note that the above regret-bound is established without assuming strong convexity of $a\mapsto -{\mathbb{E}}_{m}h\left(a,\omega \right).$ Also no Lipschitz continuity bound of the gradient is assumed.

2.4.2. Bregman learning algorithms

Algorithm 2. The Bregman learning pseudocode is as follows:

1: procedure Bregman $\left(a\left(0\right),\mathrm{ϵ},T,g,\alpha ,\beta ,m,h\right)⊳$ The Bregman learning starting from $a\left(0\right)$ within $\left[0,T\right]$

2:       $a\leftarrow a\left(0\right)$

3:       while $\mathit{\text{regret}}>\mathrm{ϵ}$ and $t\le T$ do $⊳$ We have the answer if regret is 0

4:              Compute $a\left(t\right)$ solution of (13)

5:              Compute ${\mathit{\text{regret}}}_t$

6:       end while

7:       return $a\left(t\right),{\mathit{\text{regret}}}_t$ $⊳$ get $a\left(t\right)$ and the regret

8: end procedure

Proposition 3. Let $a\mapsto {\mathbb{E}}_{m}h\left(a,\omega \right):{\mathbb{R}}^{n+2}\to \mathbb{R}$ be a concave function that has a unique global maximizer $a^{\ast }.$ Assume that $a^{\ast }$ be a feasible action profile, i.e., $a^{\ast }\in \mathcal{A}.$ Let $\alpha$ and $\beta$ be two functions such that $\dot{\beta }\left(t\right)\le e^{\alpha \left(t\right)}.$ Consider the following Bregman learning algorithm

$$\begin{array}{ll}\frac{d}{\mathit{dt}}\left[g_a\left(a\left(t\right)+e^{-\alpha \left(t\right)}\dot{a}\left(t\right)\right)\right]=e^{\alpha \left(t\right)+\beta \left(t\right)}{\nabla }_a{\mathbb{E}}_{m}h\left(a\left(t\right),\omega \right),& \\ a\left(0\right)\in {\mathbb{R}}^{n+2},\dot{a}\left(0\right)\in {\mathbb{R}}^{n+2},& \end{array}$$E13

where $a\left(0\right)$ is the initial point of the algorithm and $g$ is a strictly convex function on $a.$ Let $a\left(t\right)$ be the solution to (13). Then the average regret within $\left[t_0,T\right],t_0>0$ is bounded above by

$${\mathit{\text{regret}}}_T\le \frac{c_0}{T-t_0}{\int }_{t_0}^T{e}^{-\beta \left(s\right)}\mathit{ds},$$E14

where $c_0≔d_g\left(a^{\ast },a\left(0\right)\right)+e^{-\alpha \left(0\right)}\dot{a}\left(0\right))+e^{\beta \left(0\right)}{\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a\left(0\right),\omega \right)\right]>0.$

Proof. Let $W\left(a,\dot{a},t,a^{\ast }\right)=d_g\left(a^{\ast },a\left(t\right)+e^{-\alpha \left(t\right)}\dot{a}\left(t\right)\right)+e^{\beta \left(t\right)}{\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a\left(t\right),\omega \right)\right].$ It is clear that $W$ is positive. Moreover, $\frac{d}{\mathit{dt}}W\left(a\left(t\right),\dot{a}\left(t\right),t,a^{\ast }\right)\le 0$ for $\dot{\beta }\le e^{\alpha }.$ Thus $W\left(a\left(t\right),\dot{a}\left(t\right),t,a^{\ast }\right)\le W\left(a\left(0\right),\dot{a}\left(0\right),0,a^{\ast }\right)=c_0.$ By integration between $\left[t_0,T\right]$ it follows

$$\frac{1}{T-t_0}{\int }_{t_0}^T{\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a\left(t\right),\omega \right)\right]\mathit{dt}\le \frac{c_0}{T-t_0}{\int }_{t_0}^T{e}^{-\beta \left(s\right)}\mathit{ds}.$$

This completes the proof.

In particular, for $\beta \left(s\right)=-s+e^s,$ one obtains an error bound to the minimum value as

$$\frac{c_0}{t}{\int }_0^t{e}^{-\beta \left(s\right)}\mathit{ds}=\frac{c_0}{t}{\int }_0^t{e}^s{e}^{-e^s}\mathit{ds}=\frac{c_0\left(\frac{1}{e}-e^{-e^t}\right)}{t},$$

and for $\beta \left(s\right)=s,$ the regret bound becomes

$$\frac{c_0}{t}{\int }_0^t{e}^{-\beta \left(s\right)}\mathit{ds}=\frac{c_0\left(1-e^{-t}\right)}{t}.$$

Figure 1 illustrates the advantage of algorithm (13) compared with the gradient flow (10). It plots the regret bound $\frac{c_0}{T-t_0}{\int }_{t_0}^T{e}^{-\beta \left(s\right)}\mathit{ds}$ for $\beta =s$ and $d_g\left(a^{\ast },a_0\right)\frac{log\frac{T}{t_0}}{T-t_0}$ with an initial gap of $c_0=25.$

Figure 1.

Global regret bound under Bregman vs. gradient. The initial gap is $c_0=25.$

The advantage of algorithms (10) and (13) is that it is not required to compute the Hessian of ${\mathbb{E}}_{m}h\left(a,\omega \right)$ as it is the case in the Newton scheme. As a corollary of Proposition 2 the regret vanishes as $T$ grows. Thus, it is a no-regret algorithm. However, Algorithm (10) may not be sufficiently fast. Algorithm (13) provides a higher order convergence rate by carefully designing $\left(\alpha ,\beta \right).$ The average regret decays very quickly to zero [7]. However, it may generate an oscillatory trajectory with a big magnitude. The next subsection presents risk-aware algorithms that reduce the oscillatory phase of the trajectory.

2.4.3. Risk-aware Bregman learning algorithm

In order to reduce the oscillatory phase, we introduce a risk-aware Bregman learning algorithm [7] which is a speed-up-and-average version of (13) called mean dynamics $\overline{m}$ of $a$ given by

$$\begin{array}{ll}\stackrel{⃛}{\overline{m}}=-\frac{3}{t}\ddot{\overline{m}}-\left(e^{\alpha }-\dot{\alpha }\right)\left(\ddot{\overline{m}}+\frac{2}{t}\dot{\overline{m}}\right)& \\ +\frac{e^{2\alpha +\beta }}{t}{g}_{\overline{m}\overline{m}}^{-1}\left(\overline{m}+\left[t+2e^{-\alpha }\right]\dot{\overline{m}}+{\mathit{te}}^{-\alpha }\ddot{\overline{m}}\right)\mathbb{E}{h}_{\overline{m}}\left(t\dot{\overline{m}}+\overline{m},\omega \right),& \end{array}$$E15

with starting vector $\overline{m}\left(0\right)=a\left(0\right),\dot{\overline{m}}\left(0\right),\ddot{\overline{m}}\left(0\right).$

Algorithm 3. The risk-aware Bregman learning pseudocode is as follows:

1: procedure Risk-aware Bregman $\left(\overline{m}\left(0\right),\mathrm{ϵ},T,g,\alpha ,\beta ,m,h\right)⊳$ The risk-aware Bregman learning starting from $\overline{m}\left(0\right)$ within $\left[0,T\right]$

2:       $\overline{m}\leftarrow \overline{m}\left(0\right)=a\left(0\right),\dot{\overline{m}}\left(0\right),\ddot{\overline{m}}\left(0\right)$

3:       while $\mathit{\text{regret}}>\mathrm{ϵ}$ and $t\le T$ do $⊳$ We have the answer if regret is 0

4:             Compute $\overline{m}\left(t\right)$ solution of (15)

5:             Compute regret

6:       end while

7:       return $\overline{m}\left(t\right),{\mathit{\text{regret}}}_t$ $⊳$ get $\overline{m}\left(t\right)$ and the regret

8: end procedure

Proposition 4. The time-average trajectory of the learning algorithm (13) generates the mean dynamics (15).

Proof. We use the average relation $\overline{m}\left(t\right)=\frac{1}{t}{\int }_0^{t}a\left(s\right)\mathit{ds}$ where $a$ solves Eq. (13). From the definition of $\overline{m},$ and by Hopital’s rule, $\overline{m}\left(0\right)=a\left(0\right).$ Moreover, $\overline{m}\left(t\right)$ and $a\left(t\right)$ share the following equations:

$$\begin{array}{ll}a\left(t\right)=\overline{m}\left(t\right)+t\dot{\overline{m}}\left(t\right),& \\ \dot{a}\left(t\right)=2\dot{\overline{m}}\left(t\right)+t\ddot{\overline{m}}\left(t\right),& \\ \ddot{a}\left(t\right)=3\ddot{\overline{m}}\left(t\right)+t\stackrel{⃛}{\overline{m}}\left(t\right).& \end{array}$$E16

Substituting these values in Eq. (13) yields the mean dynamics (15). This completes the proof.

The risk-aware Bregman dynamics (15) generates a less oscillatory trajectory due to its averaging nature. The next result provides an accuracy bound for (15).

Proposition 5. The risk-aware Bregman dynamics (15) satisfies

$$0\le {\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(\overline{m}\left(t\right),\omega \right)\right]\le \frac{c_0}{t}{\int }_0^t{e}^{-\beta \left(s\right)}\mathit{ds}.$$

Proof. Let $\overline{m}\left(t\right)=\frac{1}{t}{\int }_0^{t}a\left(s\right)\mathit{ds}.$ Then, $\overline{m}\left(t\right)={\int }_{\mathbb{R}}a\left(s\right)\left(\frac{1}{t}1{\mathrm{l}}_{\left[0,t\right]}\left(s\right)\right)\mathit{ds}.$ Thus, $\overline{m}\left(t\right)={\mathbb{E}}_{\mu \left(t\right)}a$ where $\mu \left(t\right)$ is the measure with density $\mathit{d\mu }\left(t\right)\left[s\right]=\frac{1}{t}1{\mathrm{l}}_{\left[0,t\right]}\left(\mathit{ds}\right).$ By convexity of $-{\mathbb{E}}_{m}h\left(a,\omega \right)$ we apply the Jensen’s inequality:

$${\mathbb{E}}_{m}h\left(\frac{1}{t}{\int }_0^{t}a\left(s\right)\mathit{ds},\omega \right)={\mathbb{E}}_{m}h\left(\overline{m}\left(t\right),\omega \right)={\mathbb{E}}_{m}h\left({\mathbb{E}}_{\mu \left(t\right)}a,\omega \right)$$

$$\ge {\mathbb{E}}_{\mu \left(t\right)}{\mathbb{E}}_{m}h\left(a,\omega \right)=\frac{1}{t}{\int }_0^t{\mathbb{E}}_{m}h\left(a\left(s\right),\omega \right)\mathit{ds}.$$

In view of (14) one has

$$0\le {\mathbb{E}}_{m}h\left(a^{\ast },\omega \right)-{\mathbb{E}}_{m}h\left(\frac{1}{t}{\int }_0^{t}a\left(s\right)\mathit{ds},\omega \right)$$

$$\le \frac{1}{t}{\int }_0^t\left[{\mathbb{E}}_{m}h\left(a^{\ast },\omega \right)-{\mathbb{E}}_{m}h\left(a\left(s\right),\omega \right)\right]\mathit{ds}$$

$$\le c_0\frac{1}{t}{\int }_0^t{e}^{-\beta \left(s\right)}\mathit{ds},$$

$$0\le {\mathbb{E}}_{m}h\left(a^{\ast },\omega \right)-{\mathbb{E}}_{m}h\left(\overline{m}\left(t\right),\omega \right)\le \frac{c_0}{t}{\int }_0^t{e}^{-\beta \left(s\right)}\mathit{ds}.$$

This completes the proof.

Definition 5. (Convergence time). Let $\delta >0$ and $a\left(t\right)$ be the trajectory generated by Bregman algorithm starting from $a_0$ at time $t_0.$ The convergence time to be within a ball $B\left({\mathbb{E}}_{m}h\left(a^{\ast },\omega \right),\delta \right)$ of radius $\delta >0$ from the center $r\left(a^{\ast }\right)$ is given by

$$T_{\delta }=inf\left\{t|{\mathbb{E}}_m\left[h\left(a^{\ast },\omega \right)-h\left(a\left(t\right),\omega \right)\right]\le \delta ,t>t_0\right\}.$$

Proposition 6. Under the assumptions above, the error generated by the algorithm is at most (14) which means that it takes at most $T_{\delta }={\beta }^{-1}\left[log\frac{c_0}{\delta }\right]$ time units to the algorithm to be within a ball $B\left(r\left(a^{\ast }\right),\delta \right)$ of radius $\delta >0$ from the center ${\mathbb{E}}_{m}h\left(a^{\ast },\omega \right)$.

Proof. The proof is immediate. For $\delta >0$ the average regret bound of Proposition 5,

$${\mathit{\text{regret}}}_T\le \frac{c_0}{T-t_0}{\int }_{t_0}^T{e}^{-\beta \left(s\right)}\mathit{ds}\le \delta ,$$E17

provides the announced convergence time bound. This completes the proof.

See Table 1 for detailed parametric functions on the bound $T_{\delta }.$

Example 4. Let $f\left(y\right)=ylogy$ defined on ${\mathbb{R}}_{+}^{\ast }.$ Then, $f\left(1\right)=0,$ and derivatives of $f$ are $f^{\prime }\left(y\right)=1+logy,{f^{\prime }}^{\prime }\left(y\right)=\frac{1}{y}>0.$ The Legendre-Fenchel transform of $f$ is $f^{\ast }\left(\xi \right)=y^{\ast }=e^{\xi -1}.$ Let $\left(a_1,a_2\right)\mapsto g\left(a\right)=\parallel a{\parallel }_2^2,$ and $\left(a_1,a_2,\omega \right)\mapsto r\left(a_1,a_2,\omega \right)=-\left(1+{\sum }_{k=1}^2{\omega }_{\mathit{ik}}^2{a}_k^2\right).$ The coefficient $\omega$ distribution is unknown but a sampled empirical measure $m$ is considered to be similar to uniform distribution in $\left(0,1\right]$ with ${10}^4$ samples. We illustrate the quick convergence rate of the algorithm in a basic example and plot in Figure 2 the trajectories under standard gradient, Bregman dynamics and risk-aware Bregman dynamics (15). In particular, we observe that risk-aware Bregman dynamics (15) provides very quickly a satisfactory value. In this particular setup, we observe that the accuracy of the risk-aware Bregman algorithm (15) at $t=0.5$ will need four times ($t=2$) less than the standard Bregman algorithm to reach a similar level of error. It takes $40$ times more $\left(t=20\right)$ than the gradient ascent to reach that level. Also, we observe that the risk-aware Bregman algorithm is less oscillatory and the amplitude decays very fast compared to the risk-neutral algorithm.

Convergence	Error bound	Time-to-reach $T_{\delta }$
Triple exponential	$e^{-e^{e^t}}{c}_0$	$log\left[log\left(log\frac{c_0}{\delta }\right)\right]$
$\alpha \left(t\right)=t+e^t,\beta \left(t\right)=e^{e^t}$
Double exponential rate	$e^{-e^t}{c}_0$	$log\left(log\frac{c_0}{\delta }\right)$
$\alpha \left(t\right)=t,\beta \left(t\right)=e^t$
Exponential rate	$e^{-t}{c}_0$	$log\frac{c_0}{\delta }$
$\alpha \left(t\right)=0,\beta \left(t\right)=t$
Polynomial order $k$	$\frac{c_0}{t^k}$	$\frac{c_0^{1/k}}{{\delta }^{1/k}}$
$\alpha \left(t\right)=logk-logt,\beta \left(t\right)=klogt$

Table 1.

Convergence rate under different set of functions.

Figure 2.

Gradient ascent vs. risk-aware Bregman dynamics for $r=-\left(1+{\sum }_{k=1}^2{\omega }_k^2{a}_k^2\right).$

3. Constrained distributionally robust optimization

In the constrained case i.e., when $\mathcal{A}$ is a strict subset of ${\mathbb{R}}^{n+2},$ algorithms (10) and (13) present some drawbacks: The trajectory $a\left(t\right)$ may not be feasible, i.e., $a\left(t\right)\notin \mathcal{A}\times {\mathbb{R}}_{+}\times \mathbb{R}$ even when it starts in $\mathcal{A}.$ In order to design feasible trajectories, projected gradient has been widely studied in the literature. However, a projection into $\mathcal{A}$ at each time $t$ involves additional optimization problems and the computation of the projected gradient adds extra complexity to the algorithm. We restrict our attention to the following constraints:

$$\mathcal{A}=\left\{a\in {\mathbb{R}}^n|a_l\in \left[{\underset{¯}{a}}_l,{\overline{a}}_l\right],l\in \left\{1,\dots ,n\right\},\sum _{l=1}^n{c}_l{a}_l\le b\right\}.$$

We impose the following feasibility condition: ${\underset{¯}{a}}_l<{\overline{a}}_l,l\in \left\{1,\dots ,n\right\},c_l>0,$ ${\sum }_{l=1}^n{c}_l{\underset{¯}{a}}_l<b.$ Under this setting, the constraint set $\mathcal{A}$ is non-empty, convex and compact.

We propose a method to compute a constrained solution that has a full support (whenever it exists). We do not use the projection operator. Indeed we transform the domain $\left[{\underset{¯}{a}}_l,{\overline{a}}_l\right]=\xi \left(\left[0,1\right]\right)$ where $\xi \left(x_l\right)={\overline{a}}_l{x}_l+{\underset{¯}{a}}_l\left(1-x_l\right)=a_l.$ $\xi$ is a one-to-one mapping and

$$x_l={\xi }^{-1}\left(a_l\right)=\frac{a_l-{\underset{¯}{a}}_l}{{\overline{a}}_l-{\underset{¯}{a}}_l}\in \left[0,1\right].$$

$$\sum _{l=1}^n{c}_l\left({\overline{a}}_l-{\underset{¯}{a}}_l\right)x_l\le b-\sum _{l=1}^n{c}_l{\underset{¯}{a}}_l≕\widehat{b}.$$

The algorithm (18)

$$\left\{\begin{array}{l}\dot{y}={\left[{\nabla }^{2}g\right]}^{-1}{\nabla }_a{\mathbb{E}}_{m}h\left(a,\omega \right)≕\widehat{f}\left(a\right),\\ a_l≔{\overline{a}}_l{x}_l+{\underset{¯}{a}}_l\left(1-x_l\right),\\ x_l=min\left(1,\frac{e^{y_l}}{{\sum }_{k=1}^n{e}^{y_k}}\frac{\widehat{b}}{\left[c_l\left({\overline{a}}_l-{\underset{¯}{a}}_l\right)\right]}\right),\\ l\in \left\{1,\dots ,n\right\},\end{array}\right.$$E18

generates a trajectory $a\left(t\right)$ that satisfies the constraint.

Algorithm 4. The constrained learning pseudocode is as follows:

1: procedure Constrained gradient $\left(a\left(0\right),\mathrm{ϵ},T,g,m,h\right)⊳$ The constrained learning algorithm starting from $a\left(0\right)$ within $\left[0,T\right]$

2:       $a\leftarrow a\left(0\right)$

3:       while $\mathit{\text{regret}}>\mathrm{ϵ}$ and $t\le T$ do $⊳$ We have the answer if regret is 0

4:             Compute $a\left(t\right)$ solution of (18)

5:             Compute $\mathit{\text{regret}}$

6:       end while

7:       return $a\left(t\right),{\mathit{\text{regret}}}_t$ $⊳$ get a(t) and the regret

8: end procedure

Proposition 7. If $\widehat{b}\le {min}_l{c}_l\left({\overline{a}}_l-{\underset{¯}{a}}_l\right)$ then Algorithm (18) reduces to

$$\left\{\begin{array}{ll}{a}_l≔{\overline{a}}_l{x}_l+{\underset{¯}{a}}_l\left(1-x_l\right),& \\ {\dot{x}}_l=x_l\left[⟨e_l,\widehat{f}\left(a\right)⟩-\frac{1}{\widehat{b}}\sum _l⟨e_l,\widehat{f}\left(a\right)⟩x_l\left[c_l\left({\overline{a}}_l-{\underset{¯}{a}}_l\right)\right]\right],& \\ l\in \left\{1,\dots ,n\right\}& \end{array}\right.$$E19

Proof. It suffices to check that for $\widehat{b}\le {min}_l{c}_l\left({\overline{a}}_l-{\underset{¯}{a}}_l\right),$ the vector $z$ defined by $z_l=\frac{e^{y_l}}{{\sum }_{k=1}^n{e}^{y_k}}$ solves the replicator equation,

$${\dot{z}}_l=z_l\left[{\dot{y}}_l-⟨z,\dot{y}⟩\right].$$

Thus, $x_l=\frac{e^{y_l}}{{\sum }_{k=1}^n{e}^{y_k}}\frac{\widehat{b}}{\left[c_l\left({\overline{a}}_l-{\underset{¯}{a}}_l\right)\right]}$ solves ${\dot{x}}_l=x_l\left[⟨e_l,\widehat{f}\left(a\right)⟩-\frac{1}{\widehat{b}}{\sum }_l⟨e_l,\widehat{f}\left(a\right)⟩x_l\left[c_l\left({\overline{a}}_l-{\underset{¯}{a}}_l\right)\right]\right].$ This completes the proof. 

Note that the dynamics of $x$ in Eq. (19) is a constrained replicator dynamics [8] which is widely used in evolutionary game dynamics. This observation establishes a relationship between optimization and game dynamics and explains that the replicator dynamics is the gradient flow of the (expected payoff) under simplex constraint.

The next example illustrates a constrained distributionally robust optimization in wireless communication networks.

Example 5 (Wireless communication). Consider a power allocation problem over $n$ medium access channels. The signal-to-interference-plus-noise ratio (SINR) is

$${\mathit{\text{SINR}}}_l=\frac{\frac{a_l{\left|{\omega }_{\mathit{ll}}\right|}^2}{{\left(d^2\left(s_r\left(l\right),s_t\left(l\right)\right)+{\epsilon }^2\right)}^{\frac{o}{2}}}}{N_0\left(s_r\left(l\right)\right)+I_l\left(s_r\left(l\right)\right)},$$

where

• $N_0>0$ is the background noise.

• The interference on channel $l$ is denoted $I_l\ge 0.$ One typical model for $I_l$ is

$I_l=\sum _{k\ne l}\frac{a_k{\left|{\omega }_{\mathit{kl}}\right|}^2}{{\left(d^2\left(s_r\left(l\right),s_t\left(k\right)\right)+{\epsilon }^2\right)}^{\frac{o}{2}}}.$

• $\mathrm{ϵ}>0$ is the height of the transmitter antenna.

• ${\omega }_{\mathit{ll}}$ is the channel state at $l.$ The channel state is unknown. Its true distribution is also unknown.

• $s_r\left(l\right)$ is the location of the receiver of $l$

• $s_t\left(l\right)$ is the location of the transmitter of $l$

• $o\in \left\{\mathrm{2,3,4}\right\}$ is the pathloss exponent.

• $a_l$ is the power allocated to channel $l.$ It is assumed to be between ${\underset{¯}{a}}_l\ge 0$ and ${\overline{a}}_l$ with $0\le {\underset{¯}{a}}_l<{\overline{a}}_l<+\infty .$ Moreover, a total power budget constraint is imposed ${\sum }_{l=1}^n{a}_l\le \overline{a}$ where $\overline{a}>{\sum }_{l=1}^n{\underset{¯}{a}}_l\ge 0.$

It is worth mentioning that the action constraint of the power allocation problem are similar to the ones analyzed in Section 3. The admissible action space is

$$\mathcal{A}≔\left\{a\in {\mathbb{R}}_{+}^n:{\underset{¯}{a}}_l\le a_l\le {\overline{a}}_l,\sum _{l=1}^n{a}_l\le \overline{a}\right\}.$$

Clearly, $\mathcal{A}$ is a non-empty convex compact set. The payoff function is the sum-rate $r\left(a,\omega \right)={\sum }_{l=1}^n{W}_{l}log\left(1+{\mathit{\text{SINR}}}_l\right)$ where $W_l>0.$ The mapping $\left(a,\omega \right)\mapsto r\left(a,\omega \right)$ is continuously differentiable.

• Robust optimization is too conservative: Part of the robust optimization problem [9, 7] consists of choosing the channel gain ${\left|{\omega }_{\mathit{ll}}\right|}^2\in \left[0,{\overline{\omega }}_{\mathit{ll}}\right]$ were the bound $\overline{\omega }$ need to be carefully designed. However the worst case is achieved when the channel gain is zero: ${inf}_{\omega \in {\prod }_l\left[0,{\overline{\omega }}_{\mathit{ll}}\right]}r\left(a,\omega \right)=0.$ Hence the robust performance is zero. This is too conservative as several realizations of the channel may give better performance than zero. Another way is to re-design the bounds ${\underset{¯}{\omega }}_{\mathit{ll}}$ and ${\overline{\omega }}_{\mathit{ll}}.$ But if ${\underset{¯}{\omega }}_{\mathit{ll}}>0$ it means that very low channel gains are not allowed, which may be too optimistic. Below we use the distributional robust optimization approach which eliminates this design issue.

• Distributional robust optimization: By means of the training sequence or channel estimation method, a certain (statistical) distribution $m$ is derived. However $m$ cannot be considered as the true distribution of the channel state due to estimation error. The true distribution of $\omega$ is unknown. Based on this observation, an uncertainty set $B_{\rho }\left(m\right)$ with radius $\rho \ge 0$ is constructed for alternative distribution candidates. Note that $\rho =0$ means that $B_0\left(m\right)=\left\{m\right\}.$ The distributional robust optimization problem is ${sup}_a{inf}_{\tilde{m}\in B_{\rho }\left(m\right)}{\mathbb{E}}_{\tilde{m}}r\left(a,\omega \right).$ In presence of interference, the function $r\left(a,\omega \right)$ is not necessarily concave in $a.$ In absence of interference, the problem becomes concave.

4. Distributed optimization

This section presents distributed distributionally robust optimization problems over a direct graph. A large number of virtual agents can potentially choose a node (vertex) subject to constraint. The vector $a$ represents the population state. Since $a$ has $n$ components, the graph has $n$ vertices. The interactions between virtual agents are interpreted as possible connections of the graph. Let us suppose that the current interactions are represented by a directed graph $\mathcal{G}=\left(\mathcal{L},\mathcal{E}\right)$, where $\mathcal{E}\subseteq {\mathcal{L}}^2$ is the set of links representing the possible interaction among the proportion of agents, i.e., if $\left(l,k\right)\in \mathcal{E}$, then the component $l$ of $a$ can interact with the $k-$th component of $a$. In other words, $\left(l,k\right)\in \mathcal{E}$ means that virtual agents selecting the strategy $l\in \mathcal{L}$ could migrate to strategy $k\in \mathcal{L}.$ Moreover, $\mathrm{\Lambda }\in {\left\{0,1\right\}}^{n\times n}$ is the adjacency matrix of the graph $\mathcal{G}$, and whose entries are ${\lambda }_{\mathit{lk}}=1$, if $\left(l,k\right)\in \mathcal{E}$; and ${\lambda }_{\mathit{lk}}=0$, otherwise.

Definition 6. The distributionally robust fitness function is the marginal distributionally robust payoff function. If $a\mapsto {\mathbb{E}}_{m}h\left(a,\omega \right)$ is continuously differentiable, the distributionally robust fitness function is ${\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right).$

Definition 7. The virtual population state $a$ is an equilibrium if $a\in \mathcal{A}$ and it solves the variational inequality

$$⟨a-b,{\mathbb{E}}_m{\nabla }_{a}h\left(a,\omega \right)\ge 0,\forall b\in \mathcal{A}.$$

Proposition 8. Let the set of virtual population state $\mathcal{A}$ be non-empty convex compact and $b\mapsto {\mathbb{E}}_m\nabla h\left(b,\omega \right)$ be continuous. Then the following conditions are equivalent:

• $⟨a-b,{\mathbb{E}}_m\nabla h\left(a,\omega \right)\ge 0,\forall b\in \mathcal{A}.$

• the action $a$ satisfies $a={\mathit{\text{proj}}}_{\mathcal{A}}\left[a+\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)\right]$

Proof. Let $a$ be a feasible action that solves the variational inequality:

$$⟨a-b,{\mathbb{E}}_m\nabla h\left(a,\omega \right)\ge 0,\forall b\in \mathcal{A}.$$

Let $\eta >0.$ By multiplying both sides by $\eta ,$ we obtain

$$⟨a-b,\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)\ge 0,\forall b\in \mathcal{A}.$$

We add the term $⟨a,b-a⟩$ to both sides to obtain the following relationships:

$$\begin{array}{ll}⟨a-b,\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)⟩\ge 0& \forall b\in \mathcal{A},\\ \iff ⟨a-b,\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)⟩+⟨a-b,-a⟩\ge ⟨a,b-a⟩& \forall b\in \mathcal{A},\\ \iff ⟨b-a,-\left[a+\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)\right]⟩+⟨a-b,-a⟩\ge 0& \forall b\in \mathcal{A},\\ \iff ⟨b-a,a-\left[a+\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)\right]⟩\ge 0& \forall b\in \mathcal{A},\\ & \end{array}$$E20

Recall that the projection operator on a convex and closed set $\mathcal{A}$ is uniquely determined by

$$z\in {\mathbb{R}}^n,z^{\prime }={\mathit{\text{proj}}}_{\mathcal{A}}\left[z\right]\iff ⟨z^{\prime }-z,b-z^{\prime }⟩\ge 0,\forall b\in \mathcal{A}.$$

Thus

$$\begin{array}{ll}⟨b-a,a-\left[a+\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)\right]⟩\ge 0,\forall b\in \mathcal{A}& \\ \iff a={\mathit{\text{proj}}}_{\mathcal{A}}\left[a+\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)\right].& \end{array}$$E21

This completes the proof.

As a consequence we can derive the following existence result.

Proposition 9. Let the set of virtual population states $\mathcal{A}$ be a non-empty convex compact and the mapping $b\mapsto {\mathbb{E}}_m\nabla h\left(b,\omega \right)$ be continuous. Then, there exists at least one equilibrium in $\mathcal{A}.$

Proof. A direct application of the Brouwer-Schauder’s fixed-point theorem which states that if $\varphi :\mathcal{A}\to \mathcal{A}$ is continuous and $\mathcal{A}$ non-empty convex compact then $\varphi$ has at least one fixed-point in $\mathcal{A}.$ Here we choose $\varphi \left(a\right)={\mathit{\text{proj}}}_{\mathcal{A}}\left[a+\eta {\mathbb{E}}_m\nabla h\left(a,\omega \right)\right].$ Clearly $\varphi \left(\mathcal{A}\right)\subseteq \mathcal{A}$ and $\varphi$ is continuous on $\mathcal{A}$ as the mapping $b\mapsto {\mathbb{E}}_m\nabla h\left(b,\omega \right)$ and the projection operator $b\mapsto {\mathit{\text{proj}}}_{\mathcal{A}}\left[b\right]$ are both continuous. Then the announced result follows. This completes the proof.

Note that we do not need sophisticated set-valued fixed-point theory to obtain this result.

Definition 8. The virtual population state $a$ is evolutionarily stable if $a\in \mathcal{A}$ and for any alternative deviant state $b\ne a$ there is an invasion barrier ${\mathrm{ϵ}}_b>0$ such that

$$⟨a-b,{\mathbb{E}}_m\nabla h\left(a+\mathrm{ϵ}\left(b-a\right),\omega \right)>0,\forall \mathrm{ϵ}\in \left(0,{\mathrm{ϵ}}_b\right).$$

The function $\mathrm{\varrho }:\mathcal{A}\times {\mathbb{R}}^n\times {\mathbb{R}}_{+}^{n\times n}\to {\mathbb{R}}^{n\times n}$ is the revision protocol, which describes how virtual agents are making decisions. The revision protocol $\mathrm{\varrho }$ takes a population state $a,$ the corresponding fitness $\nabla {\mathbb{E}}_{m}h$, the adjacency matrix $\mathrm{\Lambda }$ and returns a matrix. Therefore, let ${\mathrm{\varrho }}_{\mathit{lk}}\left(a,h,\mathrm{\Lambda }\right)$ be the switching rate from the $l^{\mathrm{th}}$ to $k^{\mathrm{th}}$ component. Then, the virtual agents selecting the strategy $l\in \mathcal{L}$ have incentives to migrate to the strategy $l\in \mathcal{L}$ only if ${\mathrm{\varrho }}_{\mathit{lk}}\left(a,h,\mathrm{\Lambda }\right)>0$, and it is also possible to design switch rates depending on the topology describing the migration constraints, i.e., ${\lambda }_{\mathit{lk}}=0\Rightarrow {\mathrm{\varrho }}_{\mathit{lk}}\left(a,h,\mathrm{\Lambda }\right)=0.$ The distributed distributionally robust optimization consists to perform the optimization problem above over the distributed network that is subject to communication restriction. We construct a distributed distributionally robust game dynamics to perform such a task. The distributed distributionally robust evolutionary game dynamics emerge from the combination of the (robust) fitness $h$ and the constrained switching rates $\mathrm{\varrho }.$ The evolution of the portion $a_l$ is given by the distributed distributional robust mean dynamics

$${\dot{a}}_l=\sum _{k\in \mathcal{L}}{a}_k{\mathrm{\varrho }}_{\mathit{kl}}\left(a,h,\mathrm{\Lambda }\right)-a_l\sum _{k\in \mathcal{L}}{\mathrm{\varrho }}_{\mathit{lk}}\left(a,h,\mathrm{\Lambda }\right),l\in \mathcal{L},$$E22

Since the distributionally robust function $h$ is obtained after the transformation from payoff function $r$ by means of triality theory, the dynamics (22) is seeking for distributed distributionally robust solution.

Algorithm 5. The distributed distributional robust mean dynamics pseudocode is as follows:

1: procedure Population-inspired algorithm $\left(a\left(0\right),\mathrm{ϵ},T,\mathrm{\varrho },g,m,h,\mathrm{\Lambda }\right)⊳$ The population-inspired learning starting from $a\left(0\right)$ within $\left[0,T\right]$

2:       $a\leftarrow a\left(0\right)$

3:       while $\mathit{\text{regret}}>\mathrm{ϵ}$ and $t\le T$ do $⊳$ We have the answer if regret is 0

4:             Compute $a\left(t\right)$ solution of (22)

5:             Compute ${\mathit{\text{regret}}}_t$

6:       end while

7:       return $a\left(t\right),{\mathit{\text{regret}}}_t$ $⊳$ get $a\left(t\right)$ and the regret

8: end procedure

The next example establishes evolutionarily stable state, equilibria and rest-point of the dynamics (22) by designing $\mathrm{\varrho }.$

Example 6. Let us consider a power system that is composed of $10$ generators, i.e., let $\mathcal{L}=\left\{1,\dots ,10\right\}$. Let $a_l\in {\mathbb{R}}_{+}$ be the power generated by the generator $l\in \mathcal{L}$. Each power generation should satisfy the physical and/or operation constraints $a_l\in \left[{\underset{¯}{a}}_l,{\overline{a}}_l\right]$, for all $l\in \mathcal{L}$. It is desired to satisfy the power demand given by $d\in \mathbb{R}$, i.e., it is necessary to guarantee that ${\sum }_{l\in \mathcal{L}}{a}_l=d,$ i.e., the supply meets the demand. The objective is to minimize the generation quadratic costs for all the generators, i.e.,

$$\mathit{Maximize}r\left(a,\omega \right)=\sum _{l\in \mathcal{L}}{r}_l\left(a_l\right)=-\sum _{l\in \mathcal{L}}\left(c_{0l}+c_{1l}{a}_l+c_{2l}{a}_l^2\right),$$

$$\mathrm{s}.\mathrm{t}.\sum _{l\in \mathcal{L}}{a}_l=d,{\underset{¯}{a}}_l\le a_l\le {\overline{a}}_l,l\in \mathcal{L},$$

where $r:{\mathbb{R}}^n\to \mathbb{R}$ is concave, and the parameters are possibly uncertain and selected as $c_{0l}=25+6l$, $c_{1l}=15+4l+{\omega }_{1l}$, $c_{2l}=5+l+{\omega }_{2l}$, and $d=20+{\omega }_{3l}$. Therefore, the fitness functions for the corresponding full potential game are given by $f_l\left(a\right)=-2a_l{c}_{2l}-c_{1l}$, for all $l\in \mathcal{L}$, and action space is given by

$$\mathcal{A}=\left\{a\in {\mathbb{R}}_{+}^n:\sum _{l\in \mathcal{L}}{a}_l=d,a_l\in \left[{\underset{¯}{a}}_l,{\overline{a}}_l\right]\right\}.$$

The distributed revision protocol is set to

$${\mathrm{\varrho }}_{\mathit{lk}}\left(a,h,\mathrm{\Lambda }\right)=\frac{{\lambda }_{\mathit{lk}}}{a_l}max\left(0,{\overline{a}}_k-a_k\right)max\left(0,a_l-{\underset{¯}{a}}_l\right)max\left(0,{\mathbb{E}}_m\left(h_k-h_l\right)\right),$$

for $a_l\ne 0.$ We evaluate four different scenarios, i.e.,

1. $\underset{¯}{a}=0_n$ and $\overline{a}=d1{\mathrm{l}}_n$,

2. ${\underset{¯}{a}}_l=0$, for all $l\in \mathcal{L}\backslash \left\{9,10\right\}$, ${\underset{¯}{a}}_9=1.1$, and ${\underset{¯}{a}}_{10}=1$; and ${\overline{a}}_l=d$, for all $l\in \mathcal{L}\backslash \left\{1,2\right\}$, ${\overline{a}}_1=3$, and ${\overline{a}}_2=2.5$,

3. Case 1 constraints and with interaction restricted to the cycle graph $\mathcal{G}=\left(\mathcal{L},\mathcal{E}\right)$ with set of links $\mathcal{E}=\left\{{\cup }_{l\in \mathcal{L}\backslash \left\{n\right\}}\left(l,l+1\right)\right\}\cup \left\{\left(n,1\right)\right\}$,

4. Case 2 constraints and with interaction restricted as in Case 3.

Figure 3 presents the evolution of the generated power, the fitness functions corresponding to the marginal costs and the total cost. For the first scenario, the evolutionary game dynamics converge to a standard evolutionarily stable state in which $\widehat{f}\left(a^{\star }\right)=c{1}_n$. In contrast, for the second scenario, the dynamics converge to a constrained evolutionarily stable state.

Figure 3.

Economic power dispatch. Evolution of the population states (generated power), fitness functions $\widehat{f}\left(a\right)=\nabla \mathbb{E}h\left(a,\omega \right)$, and the costs $-\mathbb{E}r\left(a,\omega \right)$. Figures (a)-(c) for case 1, (d)-(f) for case 2, (g)-(i) for case 3, and (j)-(l) for case 4.

4.1. Extension to multiple decision-makers

Consider a constrained game $\mathcal{G}$ in strategic-form given by

• $\mathcal{P}=\left\{1,\dots ,P\right\}$ is the set of players. The cardinality of $\mathcal{P}$ is $P\ge 2.$

• Player $p$ has a decision space $A_p\subset {\mathbb{R}}^{n_p},n_p\ge 1.$ Players are coupled through their actions and their payoffs. The set of all feasible action profiles is $\mathcal{A}\subset {\mathbb{R}}^n,$ with $n={\sum }_{p\in \mathcal{P}}{n}_p.$ Player $p$ can choose an action $a_p$ in the set ${\mathcal{A}}_p\left(a_{-p}\right)=\left\{a_p\in A_p:\left(a_p,a_{-p}\right)\in \mathcal{A}\right\}.$

• Player $p$ has a payoff function $r_p:\mathcal{A}\to \mathbb{R}.$

We restrict our attention to the following constraints:

$$A_p=\left\{a_p\in {\mathbb{R}}^{n_p}|a_{\mathit{pl}}\in \left[{\underset{¯}{a}}_{\mathit{pl}},{\overline{a}}_{\mathit{pl}}\right],l\in \left\{1,\dots ,n_p\right\},\sum _{l=1}^{n_p}{c}_{\mathit{pl}}{a}_{\mathit{pl}}\le b_p\right\}$$

The coupled constraint is

$$\mathcal{A}=\left\{a\in \prod _p{A}_p,\sum _{p\in \mathcal{P}}⟨{\overline{c}}_p,a_p⟩\le \overline{b}\right\}.$$

Feasibility condition: If ${\underset{¯}{a}}_{\mathit{pl}}<{\overline{a}}_{\mathit{pl}},l\in \left\{1,\dots ,n_p\right\},c_{\mathit{pl}}>0,$ ${\sum }_{l=1}^{n_p}{c}_{\mathit{pl}}{\underset{¯}{a}}_{\mathit{pl}}<b_p,$ ${\overline{c}}_p\in {\mathbb{R}}_{>0}^{n_p}$ and ${\sum }_{p\in \mathcal{P}}⟨{\overline{c}}_p,{\underset{¯}{a}}_p⟩<\overline{b},$ the constraint set $\mathcal{A}$ is non-empty, convex and compact.

We propose a method to compute a constrained equilibrium that has a full support (whenever it exists). We do not use the projection operator. Indeed we transform the domain $\left[{\underset{¯}{a}}_{\mathit{pl}},{\overline{a}}_{\mathit{pl}}\right]=\xi \left(\left[0,1\right]\right)$ where $\xi \left(x_{\mathit{pl}}\right)={\overline{a}}_{\mathit{pl}}{x}_{\mathit{pl}}+{\underset{¯}{a}}_{\mathit{pl}}\left(1-x_{\mathit{pl}}\right)=a_{\mathit{pl}}.$ $\xi$ is a one-to-one mapping and

$$x_{\mathit{pl}}={\xi }^{-1}\left(a_{\mathit{pl}}\right)=\frac{a_{\mathit{pl}}-{\underset{¯}{a}}_{\mathit{pl}}}{{\overline{a}}_{\mathit{pl}}-{\underset{¯}{a}}_{\mathit{pl}}}\in \left[0,1\right].$$

$$\sum _{l=1}^{n_p}{c}_{\mathit{pl}}\left({\overline{a}}_{\mathit{pl}}-{\underset{¯}{a}}_{\mathit{pl}}\right)x_{\mathit{pl}}\le b_p-\sum _{l=1}^{n_p}{c}_{\mathit{pl}}{\underset{¯}{a}}_{\mathit{pl}}≕{\widehat{b}}_p.$$

The learning algorithm (23) is

$$\left\{\begin{array}{l}{\dot{y}}_p={\left[{\nabla }_p^{2}g\right]}^{-1}{\nabla }_{a_p}{r}_p\left(a,\omega \right),\\ a_{\mathit{pl}}≔{\overline{a}}_{\mathit{pl}}{x}_{\mathit{pl}}+{\underset{¯}{a}}_{\mathit{pl}}\left(1-x_{\mathit{pl}}\right),\\ x_{\mathit{pl}}=min\left(1,\frac{e^{y_{\mathit{pl}}}}{{\sum }_{k=1}^{n_p}{e}^{y_{\mathit{pk}}}}\frac{{\widehat{b}}_p}{\left[c_{\mathit{pl}}\left({\overline{a}}_{\mathit{pl}}-{\underset{¯}{a}}_{\mathit{pl}}\right)\right]}\right),\\ l\in \left\{1,\dots ,n_p\right\},\end{array}\right.$$E23

generates a trajectory $a_p\left(t\right)={\left(a_{\mathit{pl}}\left(t\right)\right)}_l$ that satisfies the constraint of player $p$ at any time $t.$

5. Notes

The work in [10] provides a nice intuitive introduction to robust optimization emphasizing the parallel with static optimization. Another nice treatment [11], focusing on robust empirical risk minimization problem, is designed to give calibrated confidence intervals on performance and provide optimal tradeoffs between bias and variance [12, 13]. $f$-divergence based performance evaluations are conducted in [11, 14, 15]. The connection between risk-sensitivity measures such as the exponentiated payoff and distributionally robustness can be found in [16]. Distributionally robust optimization and learning are extended to multiple strategic decision-making problems i.e., distributionally robust games in [17, 18].

Acknowledgments

We gratefully acknowledge support from U.S. Air Force Office of Scientific Research under grant number FA9550-17-1-0259.