these problems collectively as risk-based optimization (RBO). Generally, a (parametrized) RBO problem can be written as:

where x is a vector of optimization variables, θ is a vector of parameters whose values can change but are constant when solving the problem, and ξ is a random vector. The b * terms are constants and R * represent risk metrics [7], [8], [9] that measure the risk associated with the nominal objective and constraint functions f * that are random variables themselves since they depend on ξ . To evaluate the risk associated with f * , their distributions must be known. However, the underlying uncertainty distributions are never known. As such, an ambiguity set P is used to represent and reason about all distributions that have certain characteristics. This definition for (1) based on DRO includes SO and RO as special cases.

For example, if P = {P}, i.e., the distribution is known to be P, then a risk metric R is evaluated with respect to P. RBO problems are generally challenging to solve and may not be tractable. However, for some choices of ambiguity sets, risk metrics, and functions f * , (1) can be transformed into an equivalent deterministic optimization problem [9], [10], [11], [12]. In many cases, the transformations can be tedious, analogous to when transforming deterministic optimization problems into standard solver forms. Some modeling and parsing tools which have addressed certain classes of RBO problems include cvxstoc [13] and MSPPy [14] for SO, PYomo [15] with its extensions PySP for SO and ROmodel [16] and PyROS [17] for RO, and PICOS [18] and RSOME [19] for RO and DRO problems. Table I summarizes the types of problems considered by these tools.

A useful extension to these tools is the automatic generation of custom C code for embedded system applications and enhanced solve times [20]. This allows users to write optimization problems in mathematical terms and automatically generate optimized C code. Unlike the general-purpose optimization problem modeling tools mentioned earlier, there are much fewer tools that support code generation. Examples include TinyMPC [21], CVXGEN [20], and CVXPYgen [22]. TinyMPC is a recent package for generating embeddable code specifically for deterministic MPC problems. CVXGEN is an older tool for generating C code for linear and quadratic programs that is being superseded by CVXPYgen which supports a wide range of CVXPY Problem instances including

2475-1456 c ⃝ 2024 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. Authorized licensed use limited to: Univ of Texas at Dallas. Downloaded on March 17,2025 at 15:43:52 UTC from IEEE Xplore. Restrictions apply.

TABLE I PYTHON OPTIMIZATION PACKAGES COMPARISON

parameterized second-order cone programs. However, these tools do not have native support for RBO problems.

To streamline prototyping RBO problems, it is desirable to have a single tool that allows users to encode them using high-level functions and easily generate C code.

Contributions: We present cvxRiskOpt, a Python package for CVXPY to address this gap. Our package returns CVXPY Problem instances and CVXPY constraints for some RBO problems and risk constraints, making it easier to formulate RBO problems. Since cvxRiskOpt extends the CVXPY ecosystem, we can utilize CVXPYgen to generate C code as well. The main contributions of cvxRiskOpt are:

• supporting rapid prototyping with certain SO and DRO problems, including variations of chance-constrained linear programs with moment-based ambiguity sets variations of worst-case expectation problem with Wasserstein-based ambiguity sets • integrating the deterministic reformulations of the above problems into CVXPY Problem instances • utilizing CVXPYgen to generate C code. To the best of our knowledge, none of the tools that overlap with cvxRiskOpt's utility enable easily encoding RBO problems and automatic code generation, simultaneously. RSOME supports many RO and DRO problems; however, it has a few disadvantages and limitations. The modeling of ambiguity sets is based on Event-wise Ambiguity Sets [23]. While [23] provides examples of formulating several ambiguity sets as event-wise ambiguity sets (e.g., Generalized-Moment Ambiguity Sets and Wasserstein Ambiguity Sets), the formulation is complicated, which poses an additional challenge for non-experts who may simply want to experiment with some standard formulations. Compared to RSOME, PICOS provides a more user-friendly interface to a few types of RO and DRO problems that it supports. However, neither RSOME nor PICOS supports automatic code generation. While cvxRiskOpt is currently more limited in its support for RBO problems, it provides a user-friendly and beginnerfriendly interface for rapid prototyping with some widely-used risk-based problems and constraints, it extends the features of CVXPY, and it benefits from CVXPYgen for C code generation which can be used for embedded applications or faster solve times [22].

Notation: Let ⟨x, y⟩ = x T y denote the inner product, [a : b] denote set of natural numbers from a to b (inclusive), and ∥•∥ ⋆ denote the dual norm.

Our package, cvxRiskOpt, builds on top of CVXPY and provides tools to convert some convex RBO problems and constraints into CVXPY Problem classes and constraints. Upon obtaining a CVXPY Problem, its parameters can be set and the class's methods can be used to solve the problem and get an optimal solution using the CVXPY API [24]. For linear, quadratic, and second order cone Problem instances, C code can be generated using CVXPYgen.

We consider some prototypical convex RBO problems. Currently, cvxRiskOpt supports reformulating the worst-case expectation with a Wasserstein-based ambiguity set [25], into a CVXPY Problem instance (Section III) and reformulating some chance constraint linear programs (CCLPs), including the distributionally robust value-at-risk (DR-VaR) with a moment-based ambiguity set, into CVXPY constraints (Section IV).

PROBLEMS Our proposed package, cvxRiskOpt, currently supports two special classes of the general RBO problem (1). We begin with Wasserstein worst-case expectation (WCE) problems.

Consider a random variable ξ for which we have N samples ξi , i ∈ [1 : N]. Let the empirical distribution formed by these samples be PN . As the true distribution of ξ is unknown, consider a Wasserstein-based ambiguity set to be robust against, i.e., the set of all distributions in the Wasserstein ball around the empirical distribution B ϵ ( PN ) per [25, (6)]. This worst case expectation (WCE) problem [25, (10)] with respect to the Wasserstein-based ambiguity set is given by: WCE sup

where ℓ(ξ ) is a decision-independent loss function. Note that (2) may appear in (1) in the objective function (e.g., Section V-A) or as the left-hand-side of the inequality constraints (e.g., [26,Def. 3]). For a piecewise affine loss function ℓ(ξ ), (2) can be reformulated into a computationally tractable convex optimization problem [25]. In particular, here we focus on the case where ℓ(ξ ) is a piecewise affine loss function consisting of a max of affine terms. Suppose that the uncertainty set is a polytope characterized by Cξ ≤ d where C, d are a matrix and vector of appropriate dimensions. Let ℓ(ξ, x, y) = max k≤K ℓ k (ξ, x, y) with

where a k (x) is a vector that may depend linearly on x, and b k (y) = ⟨b k , y⟩ + c k is a scalar and may be linear in some decision variable y. Note that (2) is not evaluated with respect to x, y. Instead, x, y can be optimized for by another optimization problem (e.g., min x,y WCE).

From [25, (15a)], the WCE problem (2) where ℓ is a max of affine terms is equivalent to: We also implement some special cases of (3) that we also provide cvxRiskOpt classes for.

The Distributionally Robust Expectation (DR-E) is implemented as a special case of (3) where the loss function only has a single term: ℓ(ξ, x, y)

= ℓ 1 (ξ, x, y). In cvxRiskOpt, we provide the following class that inherits from WassWCEMaxAffine:

The class also supports some problem arithmetic between problems that use the same Wasserstein ball radius ϵ and same samples ξ .

Adding a WassDRExpectation and a WassWCEMax Affine problem returns an instance of WassWCEMaxAffine.

Subtraction. Subtracting WassDRExpectation problems is similar to their addition by the linearity of expectation.

Multiplication. Scalar multiplication is allowed with a WassDRExpectation problem.

Division. Scalar division of WassDRExpectation by a non-negative scalar ρ is defined as a multiplication with 1/ρ.

The Distributionally Robust Conditional Value at Risk (DR-CVaR) of a scalar random variable ξ is a special case of (3) as seen when using the formal definition of CVaR [27]:

where users provide the lists a_k= a k (x), b_k= b k (y) (here k = 2) analogously to the WassDRExpectation function.

The class also supports some problem arithmetic, particularly multiplication and division with a positive scalar.

In this section we consider two forms of the generalized linear chance constraint, using a value-at-risk (VaR) risk metric, based on [28, (6)]:

where vector a(x) is linear in a decision variable x, scalar b(y) is affine in the decision variable y and ξ 1 (z), ξ 2 are a random vector and random variable respectively and z is a decision variable. Note that the two forms are similar with the distinction being that only a or ξ 1 may depend on a decision variable, but not both since otherwise the problem would be bilinear in x, z and would no longer be convex. For ease of exposition, below we will use the common formulation P(a(x) T ξ 1 (z) + b(y) + ξ 2 ≤ 0) ≥ 1ϵ while assuming that ϕ(x, y, z) := a(x) T ξ 1 (z)+b(y)+ξ 2 is linear. Here, ϵ ∈ (0, 0.5]. Problems of the form (5) and inf

can appear as inequality constraints in (1). To aid with incorporating them, cvxRiskOpt provides functions that return deterministic linear constraints based on the reformulations in [28]. The resulting constraints can be used as constraints in a CVXPY problem. The CVXPY expression can be extracted from the CVXPY Inequality constraint (cstr.expr) to use in the objective function. Let ξ(z) = ξ 1 (z) T ξ 2 T and ã(x) = a(x) T 1 T so that

where is independent of z. We have:

We implement the reformulation from [28, Sec. 2] where it is shown that CCLP (5) ⇐⇒ κσ (x) + φ(x, y, z) ≤ 0 ( 6 ) where κ is a variable that depends on the distribution. Particularly, for Gaussian distributions, κ = -1 N (1-ϵ) where

N is the inverse cumulative distribution function of the

Authorized licensed use limited to: Univ of Texas at Dallas. Downloaded on March 17,2025 at 15:43:52 UTC from IEEE Xplore. Restrictions apply. standard Gaussian random variable [28, Sec. 2.1]. To generate the equivalent deterministic constraint, we define the function: cclp_gauss(eps, a, b, xi1_hat, xi2_hat, gam11, gam12, gam22) where xi1_hat = ξ1 , xi2_hat = ξ2 , gam11 = 11 , gam12 = γ 12 , gam22 = γ 22 .

We implement the DR-CCLP reformulation from [28,Sec. 3] where it is shown that:

for some choices of P, the ambiguity set, and κ is a variable that depends on the characteristics of P.

In particular, consider a moment-based ambiguity set with known mean and covariance P := ( ξ, ); i.e., the set of all distributions with a given mean and covariance. In this case,

28, Sec. 3.1]. If P only includes distributions that are centrally symmetric [28, Definition 3.1], then κ = 1 2ϵ . Both of these reformulations are implemented through: cclp_dro_mean_cov(eps, a, b, xi1_hat, xi2_hat, gam11, gam12, gam22, cent_symm)

where cent_symm = 1 indicates that P only includes centrally symmetric distributions.

Here, we provide some examples of using cvxRiskOpt with application to portfolio optimization, model predictive control (MPC), and moving horizon estimation (MHE). Additional examples, can be found under the "Examples" section in the package documentation: https://tsummerslab.github.io/cvxRiskOpt/. The package source code can be found at: https://github.com/TSummersLab/cvxRiskOpt. Numerical results were obtained using a MacBook Pro with an M3 Pro 12-core CPU and 18GB of RAM.

Consider the distributionally robust portfolio optimization problem from [25, Section 7.1] with yearly returns ξ of m assets and portfolio percentage weights vector x with x i ≥ 0,

This problem can be encoded as a worst-case expectation optimization problem where the objective is a max of affine terms as shown in [25, (27)] where i ∈ [1 : N] and k ∈ {1, 2}:

Listing 1. Way 1 -cvxRiskOpt implementation of (8).

Listing 2. Way 2 -cvxRiskOpt implementation of (8).

Listing 3. Solving the problem.

with

of appropriate dimensions), and the 1-norm so that ∥•∥ * is the ∞-norm.

The portfolio optimization problem (8) is convex and can be encoded using CVXPY. With cvxRiskOpt, (7) can be encoded in one of two ways.

The first uses WassDRExpectation and WassDRCVaR and the cvxRiskOpt problem arithmetic allowing users to encode the problem as shown in Listing 1 (notice how line 8 is similar to (7)). The second uses WassWCEMaxAffine to encode the reformulation (8) as in Listing 2. In both cases, i x i = 1 and x i ≥ 0 need to be added to the problem. This can be done by creating a CVXPY Problem with no objective and the two constraints, as in Listing 1 line 7, then adding that to the cvxRiskOpt problem (line 8). After setting up the problem, the parameters can be updated and the problem is solved using the chosen solver (e.g., CLARABEL) as shown in Listing 3.

Furthermore, the resulting Problem can easily be compiled into C code using CVXPYgen. Fig. 1 compares the compute times for solving (7) 1) using CVXPY only, 2) using the cvxRiskOpt functions through Listing 1, and 3) using CVXPYgen on the CVXPY Problem in Listing 1 line 8 then calling the generated C code using the Python wrapper. We refer the readers to the "Distributionally Robust Mean-CVaR Portfolio Optimization" example of the package documentation for implementation details of solving the problem with CVXPY only and compling C code. The reported time are found by aggregating the results over 100 solves. Using cvxRiskOpt adds minimal overhead and the solve times are very similar to the using CVXPY only. Compiling the code and running it using the python wrapper provides a modest speed-up of around 10%.

Consider a simple MPC-based regulator for a 1D system under process noise. The dynamics are given by: x t+1 = Ax t + Bu t + Fw t with A = 0.9, B = 0.1, F = 0.1 where w t is a Gaussian random variable w t ∼ N (w t , σ 2 w ) where w t is a time-varying mean value and σ 2 w is the variance. The stochastic MPC regulation problem starting at time t with a horizon H is given by: min

x t+τ , u t+τ are the state and control decision variables, x t is the current known state, and x min t+τ , x max t+τ are time varying bounds on the state x t+τ . Notice that the state is a random variable due to the noise. The last two constraints are chance constraints that require the state to remain within the chosen bounds with probability 1ϵ for some small ϵ value. The problem (9) can be reformulated into a deterministic risktightened problem as follows:

Finding the expected value and covariance expressions (E[x t+τ ], Cov(x t+τ ), E[x T t+τ Qx t+τ ]) is done using cvxRiskOpt helper functions, and the cclp_gauss function is used to encode the chance constraints. This simplifies the process of encoding the optimization problem while still obtaining a CVXPY Problem.

The solve times for the problem using the cvxRiskOpt functions and using the Python wrapper that calls the C code generated by CVXPYgen are reported in Fig. 2 aggregated over more than 850 solves. The compiled code results in significant speed-up by over an order of magnitude. Using cvxRiskOpt makes encoding the problem easier and enables changing the assumption about the noise, e.g., to a momentbased ambiguity set, as simple as changing the function call, while still utilizing the advantages of CVXPYgen.

Consider a moving horizon estimator (MHE) applied to:

where w t , v t are white noises with known covariances whose distribution is unknown. Denote by xt the state estimate. MHE considers the measures y t over a previous time horizon N and the current time. Let the MHE objective at time t be:

where x is the optimization variable representing the estimated state. The constraints include x t-N = xt-N , i.e., starting from the previous state estimate xt-N , and any additional constraints based on knowledge of the problem data. Here, we consider two distributionally robust chance constraints on the state: inf

where x max , x min are known state limits. The chance constraints can easily be formulated using cvxRiskOpt functions similarly to the MPC example.

For the double integrator with A = 1 dt 0 1 , B = 0.5dt 2 dt , C = 1 0 , and discrete time-step dt, the MHE problem solve times using the cvxRiskOpt functions and using the Python wrapper that calls the C code generated by CVXPYgen are reported in Fig. 3 aggregated over 200 solves for various N values. Similarly to the MPC example, cvxRiskOpt simply makes encoding the problem easier while integrating directly into CVXPY. While OSQP performs the worst, the OSQP compiled code performs the best with over a 10x speed-up for N = 90. Clarabel and ECOS compiled code provide around 3.8x and 2.2x speed-ups for N = 90, respectively.

We introduced cvxRiskOpt, a package built on top of CVXPY that provides a high-level interface for encoding some DRO problems and stochastic and distributionally robust chance constraints. Our solution is easy to use and reduces the tedious manual reformulations required with some of these problems. Since cvxRiskOpt results in and integrates with CVXPY Problem's, we can also utilize CVXPYgen to produce C code for embedded system applications. Future directions include supporting a larger range of RBO problems.