1 Introduction↩︎

Reinforcement Learning (RL) is a powerful machine learning framework that enables an agent to derive effective decision-making policies by actively engaging with a given environment. In recent years, Multi-Objective Reinforcement Learning (MORL) has gained significant interest because many real-world control problems inherently involve multiple, often conflicting objectives [1]–[9]. MORL extends standard RL to handle simultaneous optimization of multiple objectives.

The central objective of MORL is to identify a policy that yields an expected cumulative return vector lying on the Pareto frontier of the set of achievable return vectors. A common strategy involves optimizing a scalarized function defined over multiple objective returns [1], [3]. This framework seeks to identify a policy \(\pi\) that maximizes a scalarized value \(f(J_1(\pi),\ldots,J_K(\pi))\), where each \(J_k(\pi)\) represents the expected discounted return for the \(k\)-th objective among \(K (\geq 2)\) objectives, and \(f: \mathbb{R}^K \to \mathbb{R}\) is a non-decreasing scalarization function such that \(J_k(\pi) \geq J_k(\pi'), 1 \leq k \leq K \Rightarrow f(J(\pi)) \geq f(J(\pi'))\). Thus, \(f\) plays a key role in imposing the designer’s preference among multiple objectives [10], [11].

Although much of the MORL literature employs a linear \(f\) (that is, the weighted sum: \(\max_\pi \sum_{k=1}^K w_k J_k(\pi)\) with fixed \(\{ w_k \}_{k=1}^K\)) due to its simplicity, the weighted sum does not always accurately represent the preference of a designer, especially regarding fairness among objectives [3], [12]. For instance, imagine a traffic light system managing an intersection where several roads converge with asymmetric arrival rates. Instead of simply aiming to reduce the total sum waiting time for all vehicles across the roads, the designer could prioritize fairness by minimizing the longest individual waiting time among the roads. This helps reduce localized congestion [13] and avoid severe delays for individual drivers.

As such, fairness-driven objectives frequently arise in real-world scenarios and are addressed beyond the standard weighted sum, such as max-min optimization [14] or proportionally fair optimization [15] in MORL. Specifically, if the Pareto frontier contains an equalization point at which \(J_1(\pi) = J_2(\pi) = \cdots = J_K(\pi)\) is satisfied, the max-min solution attains this point [14]. This notion of max-min fairness is valuable in many applications [14], [16]–[20], such as mitigating bottlenecks in cloud and edge resource management systems [17], [18] and analyzing societal fairness across diverse demographic groups [20].

However, because the standard max-min MORL framework is formulated as an unconstrained optimization problem, it cannot address many real-world problems that involve constraints which must be satisfied. For example, in resource allocation, a MORL-based scheduler may seek to maximize throughput and fairness across task queues while operating under a strict power consumption constraint [21], [22]. Similarly, in traffic management systems, a controller must minimize the maximum total waiting time across lanes while maintaining low greenhouse-gas emissions to ensure regional sustainability. Incorporating constraints into the max-min MORL framework therefore substantially broadens its practical applicability.

Despite this motivation, incorporating constraints into unconstrained max-min MORL is not straightforward. On the one hand, existing max-min MORL algorithms may lack solution exactness, for example by optimizing a lower bound of the original max-min objective [23], [24] or by being susceptible to inexact gradient estimates [12]. On the other hand, existing constrained single-objective RL methods [25]–[29] do not account for multi-objective reward settings with \(K \geq 2\) and fairness across those objectives. Moreover, their analyses are not directly applicable to our setting due to the non-differentiability of the max-min objective.

In this paper, we propose a MORL framework that maximizes max-min fairness among homogeneous objectives while simultaneously incorporating additional quantities as constraints. To address the non-differentiability and nonlinearity of the max-min objective, we reformulate the problem as a convex optimization over occupancy measures [30] and derive a convex program admitted by the dual problem. This reformulation enables us to handle the max-min criterion and the constraints in a unified manner. We further provide a detailed theoretical foundation, including convergence guarantees for the proposed algorithm.

Our main contributions are as follows:

\(\bullet\) We introduce a MORL framework that integrates the max-min criterion into constrained optimization and show that fairness and constraint satisfaction can be addressed in a unified manner.

\(\bullet\) We propose an iterative algorithm for constrained max-min MORL and establish its theoretical foundations, accompanied by a formal convergence analysis. We also empirically evaluate its behavior in tabular settings.

\(\bullet\) We demonstrate the practical relevance of our method through applications to simulated building thermal control, multi-objective locomotion control, and greenhouse-gas-emission-aware traffic management. Our approach consistently achieves a better balance between max-min fairness and constraint satisfaction than the considered baselines.

2 Background↩︎

A multi-objective Markov decision process (MOMDP) is represented as \(\langle \mathcal{S}, \mathcal{A}, T, \mu_0, r, \gamma \rangle\), where \(\mathcal{S}\) and \(\mathcal{A}\) are the sets of states and actions, respectively, \(T\) represents the transition probability distribution, \(\mu_0\) specifies the initial state distribution, and \(\gamma \in [0,1)\) is the discount factor. The reward function \(r: \mathcal{S} \times \mathcal{A} \rightarrow \mathbb{R}^{K+L}, ~K \geq 1, L \geq 0\) is vector-valued with its \(k\)-th element denoted by \(r^{(k)} ~ (1 \leq k \leq K+L)\) such that \(|r^{(k)}| \leq r^{(k)}_{\text{max}}\), where \(K+L\) is the total number of objectives. For simplicity, we abbreviate \(\forall s \in \mathcal{S}, ~ \forall a \in \mathcal{A}\) and \(\forall (s,a) \in \mathcal{S} \times \mathcal{A}\) as \(\forall s, ~ \forall a\) and \(\forall (s,a)\), respectively.

At each timestep, the agent selects an action \(a\) in the current state \(s\) according to its (stationary) policy \(\pi: \mathcal{S} \to \mathcal{P}(\mathcal{A})\), where \(\mathcal{P}(\mathcal{A})\) represents the set of probability distributions in the action space \(\mathcal{A}\). The occupancy measure is defined as \(\rho(s,a) := \sum_{s'} \mu_0(s') \sum_{t=0}^\infty \gamma^t \text{Pr}( s_t = s, a_t = a |s_0 = s', \pi^\rho)\) where \(\pi^\rho\) is the corresponding stationary policy induced by \(\rho\), expressed as \(\pi^\rho(a|s) = \frac{\rho(s,a)}{\sum_{a'}\rho(s,a')}\) [30]. Then, the vector return evaluated by \(\pi^\rho\) is given by \[\begin{align} J(\pi^\rho) &:=& [J_1(\pi^\rho),\cdots,J_{K+L}(\pi^\rho)]^\top \in \mathbb{R}^{K+L} \nonumber \\ &=&\mathbb{E}_{\pi^\rho} \left[ \sum_{t=0}^\infty \gamma^t r_t \right] = \sum_{(s,a)} r(s,a) \rho(s,a). \end{align}\]

In this paper, we define \(\mathbb{R}_{+}^L := \{ u \in \mathbb{R}^L | u_l \geq 0, ~ 1 \leq l \leq L \}\) and \(\Delta^K := \{ w \in \mathbb{R}^K | \sum_{k=1}^K w_k = 1; ~ w_k \geq 0, ~ 1 \leq k \leq K \}\), i.e., the \((K-1)\)-dimensional simplex.

3 Constrained Max-Min MORL Framework↩︎

3.1 Theoretical Foundation↩︎

We consider the MORL setting where the last \(L\) of the total \(K+L\) objectives should satisfy certain constraints. For theoretical development in this section, we assume that \(\mathcal{S}\) and \(\mathcal{A}\) are finite. The problem is formulated as follows: \[\begin{align} &&\max_{\pi^\rho} ~f( J_1(\pi^\rho),\cdots,J_{K}(\pi^\rho)) + \beta \sum_{s} \mathcal{H}_\rho(s) \rho(s) \tag{1} \\ && ~~s.t.~~ J_{K+l}(\pi^\rho) \geq C^{(l)}, ~~ l =1, \cdots, L \tag{2} \end{align}\] where \(\mathcal{H}_\rho(s) := - \sum_a \pi^\rho(a|s) \log \pi^\rho(a|s)\) is the entropy of \(\pi^\rho(\cdot|s)\), \(\rho(s) := \sum_a \rho(s,a)\) is the stationary state distribution in \(\mathcal{S}\), \(\beta > 0\) is a balancing coefficient, and \(\{ C^{(l)} \}_{l=1}^L\) is a set of threshold values. We assume a mild condition that the set \(\{ C^{(l)} \}_{l=1}^L\) is chosen by the designer such that the optimization in 1 and 2 is feasible, an assumption commonly made in the constrained MDP literature [26], [31].

In this paper, we set \(f\) the minimum function, i.e., \(f( J_1(\pi^\rho),\cdots,J_{K}(\pi^\rho))=\min_{1 \leq k \leq K} J_k(\pi^\rho)\). We note that the entropy term is required to eliminate the indeterminacy of the max-min solution without any regularization [12]. The problem reduces to the unregularized formulation as \(\beta \to 0\), with the optimality gap decreasing linearly:

Proposition 1 shows that the regularized problem is a valid approximation of the unregularized criterion.

Since directly optimizing 1 and 2 with \(f=\min\) and \(J_k(\pi^\rho) = \mathbb{E}_{\pi^\rho}[\sum_{t=0}^\infty \gamma^t r^{(k)}_t]\) is non-trivial due to its non-differentiable and nonlinear structure, we address this challenge using the occupancy measure (i.e., stationary distribution [30]) formulation. The above optimization problem with \(f=\min\) can be rewritten as \[\label{eq:maxmin95obj} \max_{\rho \geq 0} \min_{1 \leq k \leq K} \bigg( \sum_{(s,a)} r^{(k)}(s,a) \rho(s,a) \bigg) + \beta \sum_{s} \mathcal{H}_\rho(s) \rho(s)\tag{3}\] \[\label{eq:maxmin95flow} \sum_{a'} \rho(s',a') = \mu_0(s') + \gamma \sum_{(s,a)} T(s' | s,a) \rho(s,a) ,~\forall s'\tag{4}\] \[\label{eq:maxmin95const95c} \sum_{(s,a)} c^{(l)}(s,a) \rho(s,a) \geq C^{(l)}, ~~ l =1, \cdots, L\tag{5}\] where 4 is the Bellman flow equation for the occupancy measure [30]. Here, we use the notation \(c^{(l)}(s,a) := r^{(K+l)}(s,a), ~l=1,\cdots,L\) to explicitly represent the dimensions associated with the constraint. For example, these quantities can be negative of some costs. Then the formulation in 3 , 4 , and 5 constitutes a convex optimization problem.

Now we derive a convex optimization problem admitted by the dual of 3 , 4 , and 5 , which forms the foundation of our algorithm. Given \((u,w)\), let \(v^{*}_{u, w}\) be the fixed point of the operator \(\mathcal{T}_{u,w}\): \[\begin{align} \label{eq:equivalent95fixed95point} [\mathcal{T}_{u,w} v](s) &= \beta \log \sum_a \exp [ \frac{1}{\beta} \{ \sum_{l=1}^L u_l c^{(l)}(s,a) \nonumber \\ &+ \sum_{k=1}^K w_k r^{(k)}(s,a) + \gamma \sum_{s'}T(s'|s,a) v(s') \} ], \forall s. \end{align}\tag{6}\]

We emphasize that the optimization problem in ?? is not the dual per se of 3 , 4 , and 5 , since the dual problem itself does not involve the fixed point \(v^{*}_{u, w}\) in the objective. Consequently, establishing the convexity of ?? requires proving the convexity of \(v^{*}_{u, w}\) with respect to \((u,w)\), as explained in Appendix 8. Throughout our analysis, we assume that Slater condition holds [32], [33], which encompasses the standard constraint feasibility assumption [26], [31].

Proposition 2 hints that \(v^{*}_{u, w}\) can be obtained via 6 and the weights \((u,w)\) can be obtained by minimizing the loss \(\mathcal{L}(u,w)\) in ?? by some method. However, solving the optimization problem ?? is non-trivial because the fixed point \(v^*_{u,w}\) does not have a closed-form expression in terms of \((u,w)\). To address this issue, we derive the key properties of \(v^*_{u,w}\). We define a policy \(\pi^*_{u,w}\) as \[\label{eq:optimal95policy} \pi^*_{u,w}(a | s) := \frac{ \exp ( \frac{1}{\beta} Q^*_{u,w}(s,a) ) }{\sum_{a'} \exp ( \frac{1}{\beta} Q^*_{u,w}(s,a') ) }\tag{7}\] where \[\begin{align} \label{eq:eq10} Q^*_{u,w}(s,a) &:= \sum_{l=1}^L u_l c^{(l)}(s,a) + \sum_{k=1}^K w_k r^{(k)}(s,a) \nonumber \\ &+ \gamma \sum_{s'}T(s'|s,a) v^{*}_{u,w}(s'). \end{align}\tag{8}\]

Then, \(\pi^*_{u,w}\) is an optimal policy for the entropy-regularized RL [34] with a scalar reward \(\sum_{l=1}^L u_l c^{(l)}(s,a) + \sum_{k=1}^K w_k r^{(k)}(s,a), \forall (s,a)\). We derive the following theorem characterizing the relationship between \(\pi^*_{u,w}\) and the gradient of \(v^{*}_{u, w}\):

The main advantage of Theorem 3 lies not only in the simplicity of its gradient expression, but also in the fact that the gradient components corresponding to the unconstrained (\(w\)) and constrained (\(u\)) parts are obtained in a unified manner and are directly linked to the value function evaluated under \(\pi^*_{u,w}\). This unified form enables us to apply gradient descent to the optimization problem using the gradient \((\nabla_{v} v^*_{u,w}(s),\nabla_{w} v^*_{u,w}(s))\) combined with value iteration.

First, consider the unconstrained part. From ?? , \(\nabla_w \mathcal{L}(u,w)\) is given by \(\mathbb{E}_{s \sim \mu_0}[ v_r^{\pi^*_{u,w}}(s) ] \in \mathbb{R}^K\). Consequently, dimensions with smaller values of \(\mathbb{E}_{s \sim \mu_0}[ v_r^{\pi^*_{u,w}}(s) ]\) undergo smaller decreases in \(w_k\) under gradient descent, resulting in larger \(w_k\) relative to other dimensions after update with projection onto \(\Delta^K\). As a result, dimensions with smaller values receive larger updated weights, thereby realizing the max-min principle. The constrained component follows a similar logic to enforce feasibility, enabling the simultaneous achievement of max-min fairness and constraint satisfaction. In addition, our framework encompasses unconstrained max-min MORL (\(L=0\) with learned \(w\)) as well as both constrained and unconstrained RL with fixed \(w\) such as max-average RL (\(w=[1/K, \cdots, 1/K] \in \Delta^K\)).

A natural approach to solving the convex optimization problem in ?? , using Theorem 3, is projected gradient descent, since the variables \((u,w)\) lie in the convex set \(\mathbb{R}_+^L \times \Delta^K\). The convergence of projected gradient descent depends on the smoothness of the objective function [32], [35]. In our setting, \(\mathcal{L}(u,w)\) satisfies the following smoothness property, which is essential for establishing the convergence guarantee in Section 3.2:

In other words, \(\nabla_{(u,w)} \mathcal{L}(u,w)\) is \(\alpha\)-Lipschitz continuous in \(\| \cdot \|_2\). For the proof of Theorem 4, we first show that \(\nabla v^*_{u,w}\) is differentiable in \((u,w) \in \mathbb{R}^{L+K}\) (equivalently, that \(v^*_{u,w}\) is twice-differentiable) in Proposition 5. Then incorporating the closed-form expression of the Hessian of \(v^*_{u,w}(s)\) for each \(s\) into the generalized mean value inequality constitutes a key ingredient in the proof.

Let \(Q^{\pi^*_{u,w}}(s,a) \in \mathbb{R}^{L+K}\) be the action-value function evaluated with the policy \(\pi^*_{u,w}\), and define \(B^{\pi^*_{u,w}}(s) := \mathbb{E}_{a \sim \pi^*_{u,w}(\cdot|s)} \bigg[ ( Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot|s)} [Q^{\pi^*_{u,w}}(s,a')] ) ( Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot|s)} [Q^{\pi^*_{u,w}}(s,a')] )^\top \bigg] \in \mathbb{R}^{(L+K) \times (L+K)}\).

This result reflects the fact that the optimization problem in ?? is a convex program, as established in Proposition 2, owing to the positive semidefiniteness of ?? from \((I_p - \gamma T^{\pi^*_{u,w}})^{-1} = \sum_{i=0}^\infty (\gamma T^{\pi^*_{u,w}})^i\) and the relation \(H[\mathcal{L}(u,w)] = \mathbb{E}_{s \sim \mu_0} [H[v^*_{u,w}(s)]]\).

Note that \(Q^{\pi^*_{u,w}}(s,a)\) in Proposition 5 is different from \(Q^*_{u,w}(s,a)\) in 8 . By definition in the entropy-regularized RL, \(Q^*_{u,w}(s,a) \in \mathbb{R}\) is the cumulative scalarized return plus the cumulative entropy sum from \(\pi^*_{u,w}\). On the other hand, \(Q^{\pi^*_{u,w}}(s,a) \in \mathbb{R}^{L+K}\) is a cumulative sum of unconstrained rewards and constrained rewards from \(\pi^*_{u,w}\) without the entropy sum. Therefore, \([u;w]^\top Q^{\pi^*_{u,w}}(s,a)\) equals to \(Q^*_{u,w}(s,a)\) minus the cumulative entropy sum of \(\pi^*_{u,w}\).

3.2 Algorithm and Convergence Analysis↩︎

Based on the foundation built in the previous section, we propose an algorithm for constrained MORL with max-min fairness. Note that we need to jointly update the weights \((u,w)\) and the value function, which approximates \(v^{*}_{u, w}\). We adopt the following update method alternating between update of the value function and the weights \((u,w)\).

First, given a weight \((u,w)\), we update the value function to realize 6 . For this, we use an action value function \(Q\), which approximates \(Q^*_{u,w}\). If we plug 8 into the right-hand side of 6 , we have \(v^{*}_{u, w}(s) = [\mathcal{T}_{u,w} v^{*}_{u, w}](s) = \beta \log \sum_a \exp \left( \frac{Q^*_{u,w}(s,a)}{\beta} \right)\) for each \(s\). Using this form of \(v_{u,w}^*(s)\), we implement applying \(\mathcal{T}_{u,w}\) as updating the Q-function with the following: \[\label{eq:sql} \begin{align} Q(s,a) \leftarrow {} & [u; w]^\top [c; r] \\ & + \gamma \sum_{s'} T(s'|s,a)\, \beta \log \sum_{a'} \exp\!\left(\frac{Q(s',a')}{\beta}\right), \end{align}\tag{10}\]

for all \((s,a)\). In practice, 10 is iterated until the maximum change in the \(Q\)-function across successive iterations is below a given threshold.

Next, we compute an estimated gradient of \(\nabla_{(u,w)}\mathcal{L}(u,w)\) at the current weight \((u,w) = (u^m,w^m)\) where \(m=1,2,\cdots\) is the iteration index. We have shown that \(\nabla_{u} v^*_{u, w}(s) = v_c^{\pi^*_{u,w}}(s), ~ \nabla_{w} v^*_{u, w}(s) = v_r^{\pi^*_{u,w}}(s)\) for each \(s,\) where we denote \(v_c^{\pi^*_{u,w}}(s) \in \mathbb{R}^L\), \(v_r^{\pi^*_{u,w}}(s) \in \mathbb{R}^K\) as the value functions evaluated with the policy \(\pi^*_{u,w}\) for constrained reward \(c\) and unconstrained reward \(r\), respectively. Note that the policy is extracted from the updated Q-function in 10 based on the form of 7 . We then update \((u,w)\) using projected gradient descent: \[\label{eq:proj95grad} (u^{m+1},w^{m+1}) = \mathcal{P}_{K,L} [ (u^m,w^m) - l_{\text{w}} \nabla_{(u,w)}\mathcal{L}(u^m,w^m) ]\tag{11}\] where \(l_{\text{w}}\) is a learning rate for \((u,w)\) and \(\mathcal{P}_{K,L}[\cdot]\) is the projection onto the \(\mathbb{R}_+^L \times \Delta^K\). We use the convex optimization method from [36] to project onto the simplex \(\Delta^K\), and apply non-negativity clipping for projection onto \(\mathbb{R}_+^L\). Note that the projection onto \(\Delta^K\) is numerically stable as it is fully deterministic and avoids randomized procedures. In addition, its complexity is \(O(K \log K)\) [36] which is lightweight due to the logarithmic term.

We iterate this process for each \(m\), and Algorithm 1 presents the pseudocode of our algorithm.

We now provide our convergence analysis under the following assumption. Let \(\rho^{\pi^*_{u,w}}(s) = \sum_{t=0}^\infty \gamma^t \text{Pr}(s_t=s | \pi^*_{u,w}, T, \mu_0)\) be the discounted state visitation frequency of policy \(\pi^*_{u,w}\) at state \(s\) [30].

Assumption There exists at least one state \(s \in \mathcal{S}\) such that \(\rho^{\pi^*_{u,w}}(s) > 0\) and the centered action-value vectors \(\bigl\{\,Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot\mid s)}[\,Q^{\pi^*_{u,w}}(s,a')\,]\;:\;a\in\mathcal{A}\bigr\}\) span \(\mathbb{R}^{K+L}\).

This condition fails only in degenerate multi-objective settings when for every state \(s\in\mathcal{S}\) with \(\rho^{\pi^*_{u,w}}(s) > 0\), the set \(S_{\text{center}}(s)\) lies entirely within an affine subspace of dimension less than \(K+L\). Under Assumption, the Hessian \(H[\mathcal{L}(u,w)]\) is positive definite. We note that the property \(\pi^*_{u,w}(a|s) > 0\) for all \((s,a)\), implied by 7 , plays an essential role in ensuring the positive definiteness of \(H[\mathcal{L}(u,w)]\), which constitutes an additional technical benefit of the entropy regularization in 3 . (See Appendix 12.1 for more details.)

Moreover, under Slater condition with strict feasibility of the \(L\) constraints, we can restrict \(u \in \mathbb{R}^L_+\) to a compact set \(\mathcal{C}\) [29], ensuring the existence of \(\lambda := \min_{(u,w) \in \mathcal{C} \times \Delta^K} \lambda_{\text{min}}(H[\mathcal{L}(u,w)])\), where \(\lambda_{\min}(\cdot)\) denotes the minimum eigenvalue of a real symmetric matrix. (See Appendix 12.2 for more details.) Furthermore, \(\lambda\) satisfies \(0 < \lambda \leq \alpha\) [35].

Theorem 6 provides a formal guarantee of convergence for Algorithm 1 under approximate Q-updates.

Theorem 6 establishes that the error decreases geometrically at rate \(O\left((1-\tfrac{\lambda}{\alpha})^m\right)\), up to \(O(\epsilon)\). For completeness, Appendix 12.4 provides an analysis of the case without Assumption, supporting the practical implementation of our algorithm beyond Assumption in Section 5.2.

4 Related Work↩︎

MORL   The dominant approach in MORL is utility-based [1], [3], where the goal is to find an optimal policy \(\pi^* = \arg \max_{\pi} f(J(\pi))\) given a non-decreasing scalarization function \(f: \mathbb{R}^K \rightarrow \mathbb{R}\). When \(f\) is linear, each non-negative weight vector defines a scalarized MDP [37], motivating work on learning a single model capable of generating policies across a continuum of preferences

[2], [9], [38]–[42]. This family of approaches is known as multi-policy MORL [1], [3].

For nonlinear scalarization functions, however, Bellman optimality no longer holds in its standard form due to the loss of linearity, making optimization substantially more difficult [1], [3]. Algorithms that directly optimize a given nonlinear scalarized objective fall under the single-policy MORL category [1], [3]. Most work in this category considers a welfare function [8], [10], [11] as the nonlinear scalarization \(f\). Note that single-policy and multi-policy methods are complementary rather than interchangeable [1], [3].

Unconstrained Max-min MORL   Max-min MORL studies the case where \(f=\min\), aiming to enforce max-min fairness. This formulation is useful in many applications such as mitigating bottlenecks in cloud and edge resource management systems [17], [18]. Several studies optimize proxy objectives related to the unconstrained max-min formulation, for example, maximizing a conservative lower bound of \(\mathbb{E}_\pi \left[ \min_{1 \leq k \leq K}( \sum_{t=0}^\infty \gamma^t {r_t^{(k)}} ) \right]\) [23], [24], or maximizing the total return while enforcing per-group performance constraints [20].

[12] proposed a tractable approach for exact unconstrained max-min MORL using Gaussian smoothing to estimate gradients. However, this approach requires maintaining multiple network copies, increasing computational overhead. Furthermore, the gradient estimates are inherently inexact, as Gaussian smoothing of a convex function yields a convex upper bound rather than the true function [43]. In contrast, our method produces direct, theoretically grounded gradient estimates and extends naturally to constrained MORL. Recently, [44] introduced an alternative unconstrained max-min MORL formulation based on a game-theoretic framework [45], [46]. However, their method does not address the incorporation of constraints.

Constrained RL   Many approaches to constrained MDPs reformulate the problem with a scalar reward (i.e., a special case of 1 and 2 with \(K=1\) and without \(f\)) into an unconstrained one by augmenting the objective with a weighted sum of constraint violations, typically via a Lagrangian formulation [26]–[29], [31], [47], [48]. The motivation for this line of work is that the Lagrangian relaxation exhibits no duality gap, even when the original problem is non-convex with respect to the policy [48]. However, these approaches do not consider the multi-objective reward setting in 1 and 2 with \(K \geq 2\). Moreover, applying them directly to our setting is non-trivial, since \(f = \min\) introduces non-differentiability in 1 .

To resolve this, we reformulate our problem as a convex program using occupancy measures and then derive another convex program admitted by the dual problem. In particular, we show that both the max-min criterion and the constraints can be satisfied by jointly updating the weights \(u\) and \(w\), a simple yet effective approach that to our knowledge has not been explored in the constrained MDP literature. Although [49] also leverages convex analysis with occupancy measures, its focus is on constrained single-objective RL with a scalar reward (i.e., \(K=1\)) in an offline setting. Unlike our work, it neither addresses max-min fairness across multiple objectives nor provides a convergence analysis.

Several recent works have incorporated constraints into MORL, but under settings that differ from our framework which explicitly integrates max-min optimization. Some approaches learn preference-vector-conditioned policies to explore trade-offs among objectives [7], [50]. While complementary to our work, it remains unclear how to select a preference vector such that the resulting policy exactly corresponds to the optimal solution under nonlinear scalarizations, such as the max-min criterion. In contrast, our approach directly solves a single-policy constrained max-min optimization problem, mirroring the distinction between single- and multi-policy methods in unconstrained MORL [1], [3]. [51] study offline constrained MORL, where policies are trained on offline data and adapted to target preferences using additional demonstrations. In contrast, our work focuses on online learning and does not rely on offline data or demonstrations. [6] improve Pareto frontier coverage in unconstrained MORL by training multiple policies via constrained optimization, rather than directly addressing constrained MORL.

5 Experiments↩︎

In this section, we present experimental validations of our theoretical analysis and algorithm. Section 5.1 examines the convergence behavior of our method in tabular settings. In Section 5.2, we further demonstrate the practical relevance of our approach through applications beyond tabular settings including building thermal control, multi-objective locomotion control, and greenhouse-gas-emission-aware traffic management.

5.1 Tabular Settings↩︎

We conducted experiments in tabular settings to evaluate the convergence of our algorithm. Constrained MOMDPs were randomly generated, after a feasibility check, within two widely used classes of structured MDPs. (See Appendix 13 for details on the experimental setup, feasibility check, and baselines.) First, bipartite state graphs partition the state space into two disjoint subsets, enforcing transitions between them at alternating time steps. This structure captures temporal dynamics in systems with role alternation or interleaving phases [52]. Second, hierarchical MDPs organize the state space into multiple levels or stages, where transitions flow sequentially from one level to the next. This reflects tasks with subgoals or temporal abstraction [53].

The optimal value for each constrained MOMDP was computed by solving 3 , 4 , and 5 with \(\beta=0\) via linear programming (LP), and performance was evaluated as the error relative to these LP-optimal values. The goal of this experiment is to verify that our algorithm, which simultaneously accounts for both max-min fairness and constraint satisfaction, achieves near-optimal performance compared to methods that lack either max-min fairness or constraint enforcement. We compare our constrained max-min algorithm against unconstrained max-min, constrained max-average which does not enforce max-min fairness, and unconstrained max-average baselines.

5pt

Table 1 shows that the constrained max-min approach achieves near-optimal values compared to the other methods. Baselines that incorporate neither max-min fairness nor constraint satisfaction exhibit the largest optimal value errors. We further analyze the effect of \(\beta\) on convergence. As shown in Table 2, \(\beta=0.03\) and \(0.01\) yield the smallest errors. The slight increase in optimal value error when \(\beta < 0.01\) may be attributed to the first term in 12 to converge more slowly as the iteration count \(m \to \text{ITER}\), the maximum number of iterations used in Algorithm 1.

5.2 Extension to Applications↩︎

In this section, we extend our algorithm to practical applications, demonstrating that our method is not limited to tabular settings. The Building (Section 5.2.1) and MoAnt (Section 5.2.2) environments have continuous state and action spaces, while the Traffic environment (Section 5.2.3) has a continuous state space and a discrete action space. To ensure stable gradient estimation of our algorithm in continuous state spaces, we parameterize a gradient network \(g_\theta(s) \in \mathbb{R}^{L+K}\) to estimate \(\nabla_u v^{*}_{u, w}(s)\) and \(\nabla_w v^{*}_{u, w}(s)\), following Theorem 3. Implementation details for deep neural network policies, including gradient estimation and our constrained max-min algorithm for practical applications, are provided in Appendix 14.1. Our sourcecode can be found in https://github.com/Giseung-Park/Constrained-Maxmin-MORL.

5.2.1 Building Thermal Control↩︎

We consider a thermal control environment for a building from SustainGym [54]. The environment models a large commercial building with \(N_{\text{zone}}=23\) thermal zones, where the system controls the zone temperatures under changing outside conditions and occupancy. At each timestep, the system observes the current temperatures of the \(N_{\text{zone}}\) zones together with additional indicators (e.g., outdoor temperature, global horizontal irradiance, ground temperature, occupancy power, carbon intensity, and electricity price). Then the system outputs a continuous action that controls the heating/cooling power of each zone.

The zones are grouped into \(K=3\) categories corresponding to the bottom, middle, and top floors of the building. The unconstrained reward is a \(K\)-dimensional vector, where each component quantifies the comfort of the corresponding group based on the proximity of the zone temperatures to their target values. The cost is defined as the total power consumption of the system, normalized by the environment. The system’s objective is to control the building so as to maximize the minimum cumulative comfort across all groups, while satisfying the energy constraint with its designed threshold value \(C_{th}=180\).

We consider five baselines: (i) a random policy that samples a continuous action at each timestep (Random), (ii) unconstrained max-average SAC (MA-SAC) [55], (iii) max-average SAC with a Lagrangian relaxation (MA-SAC-L) [26], [56], (iv) unconstrained max-min MORL algorithm with Gaussian smoothing (Max-min GS) [12], and (v) unconstrained max-min MORL algorithm based on a two-player zero-sum game (ARAM) [44]. Each of the baselines lacks either max-min fairness ((iii)), constraint handling ((iv),(v)), or both ((i), (ii)). We report the mean performance computed across five random seeds. (See Appendices 14.2 and 14.3 for the implementation of the max-min GS baseline and for details on the environment and the hyperparameters, respectively.)

Table 3 presents the cumulative cost sum and the worst-group comfort return. Compared to the Random baseline, MA-SAC increases the worst-group comfort return but still fails to satisfy the power consumption constraint, with its cost sum exceeding the threshold \(C_{th}=180\). While MA-SAC-L satisfies the power constraint, it does so at the cost of a lower worst-group comfort return compared to MA-SAC. However, our method achieves a higher worst-group comfort return relative to MA-SAC-L, while still adhering to the power constraint. We note that both Max-min GS and ARAM violate the power constraint.

Table 4 shows that removing the max-min-related \(w\) update reduces the worst-group comfort return, while ablating the constraint-related \(u\) update causes constraint violations. These results confirm that our method effectively balances max-min performance with constraint satisfaction. (For further analysis of the effect of \(\beta\) and the computational complexity, see Appendix 14.4.)

5.2.2 Multi-objective Locomotion Control↩︎

We include MoAnt-v5 environment [57], where the agent learns locomotion to maximize \(x\) and \(y\) velocities while keeping energy consumption under a threshold. We consider an asymmetric case where movement in the \(x\) direction is attenuated by friction at rate 0.3. The velocities \((v_x, v_y)\), combined with bonus terms, constitute a two-dimensional reward, while the control cost is treated as a constraint. (See Appendix 14.5 for details on hyperparameters.)

Table 5 shows that both our method and MA-SAC-L satisfy the constraints, but our method achieves superior max-min performance. In contrast, the other four algorithms severely violate the constraints, as they do not explicitly account for constraint satisfaction. Overall, our algorithm balances constraint satisfaction and max-min fairness. (See Appendix 14.6 for the result with \(C_{th}=40\).)

5.2.3 Greenhouse-gas-emission-aware Traffic Management↩︎

We extend our method to an environment with a larger objective space. We note that MORL benchmark environments, particularly those with more than four objectives, remain limited [3]. Recently, [42] proposed a traffic signal control environment with a 16-dimensional objective vector to evaluate scalability in the MORL setting. We modify this environment to incorporate a setting aimed at reducing greenhouse-gas emissions while maintaining max-min fairness.

In a simulated 16-lane four-way intersection, the agent receives a 37-dimensional continuous state encoding traffic information at each timestep. The agent manages the traffic lights by choosing the traffic signal phase as its action. The feedback signal consists of a 16-dimensional reward vector, where each component represents the negative total waiting time for a corresponding lane. Following prior work [42], we evaluate performance in an asymmetric traffic flow scenario. Vehicles enter from each direction according to predefined inflow proportions and may proceed straight or make left or right turns.

Our goal is to achieve fair traffic flow across all lanes while enforcing a constraint on total \(\isotope{CO}_2\) emissions, contributing to a more sustainable traffic control system. Since this environment operates with a discrete action space, we replaced the MA-SAC and MA-SAC-L baselines with policy-gradient-optimization (PGO) variants, namely MA-PGO and MA-CPGO, respectively. (See Appendix 14.7 for implementation details.)

As shown in Table 6, both our method and MA-CPGO satisfy the constraints; however, our approach achieves better max-min performance. In contrast, the remaining four algorithms exhibit constraint violations. Overall, our method consistently balances max-min fairness with constraint satisfaction.

6 Conclusion↩︎

We have proposed a MORL framework that integrates max-min fairness with constraint satisfaction. Our approach offers flexibility in modeling problems that satisfy fairness and operational constraints. We establish a theoretical foundation and develop an algorithm that achieves a balance between max-min fairness and constraint satisfaction, as demonstrated in our experiments. A discussion of limitations and future work is provided in Appendix 15.

Acknowledgment↩︎

We sincerely thank anonymous reviewers for their valuable feedback during the review process. This work was supported by the Institute of Information & Communications Technology Planning & Evaluation (IITP) under Grant No. RS-2024-00457882 (AI Research Hub Project) and the National Research Foundation of Korea (NRF) grant (2022K1A3A1A31093462) funded by the Korea government (MSIT), and by the Israel Ministry of Innovation, Science & Technology, under Grant No. 1001556423 and the ISF grant 2197/22.

Impact Statement↩︎

In this work, we propose a MORL algorithm that integrates max-min fairness with constraint satisfaction. First, max-min MORL plays a critical role in promoting fairness across objectives in domains such as traffic management and resource allocation. Unfair results can lead to user dissatisfaction and, in turn, degrade overall system performance, for example, by contributing to traffic congestion. Second, incorporating constraints into RL systems is essential for the responsible development of AI systems, especially given real-world limitations on resources such as electricity, power consumption, and fossil fuels.

Our work advances the goal of sustainable AI by simultaneously incorporating fairness and resource constraints into decision-making. This contrasts to traditional methods that prioritize performance alone, often overlooking concerns of equity and efficient resource use. We believe our framework has the potential to make a meaningful and positive impact on the broader AI community, not only in resource allocation but also in emerging areas such as fair and safe alignment of large language models.

7 Proof on Optimality Gap↩︎

8 Proof on Equivalent Optimization↩︎

9 Proof of Differentiability↩︎

10 Proof of Smoothness↩︎

11 Proof of Twice-Differentiability↩︎

12 Convergence Analysis↩︎

12.1 Assumption for Action-value Nondegeneracy↩︎

Assumption There exists at least one state \(s \in \mathcal{S}\) such that \(\rho^{\pi^*_{u,w}}(s) > 0\) and the centered action-value vectors \(\bigl\{\,Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot\mid s)}[\,Q^{\pi^*_{u,w}}(s,a')\,]\;:\;a\in\mathcal{A}\bigr\}\) span \(\mathbb{R}^{K+L}\).

This condition fails only in degenerate multi-objective settings when for every state \(s\in\mathcal{S}\) with \(\rho^{\pi^*_{u,w}}(s) > 0\), the set \(\bigl\{\,Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot\mid s)}[\,Q^{\pi^*_{u,w}}(s,a')\,]: a\in \mathcal{A}\bigr\}\) lies entirely within an affine subspace of dimension less than \(K+L\) (e.g., the size of an action set is smaller than the number of objectives).

Then \(B^{\pi^*_{u,w}}(s) = \mathbb{E}_{a \sim \pi^*_{u,w}(\cdot|s)} \bigg[ ( Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot|s)} [Q^{\pi^*_{u,w}}(s,a')] ) ( Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot|s)} [Q^{\pi^*_{u,w}}(s,a')] )^\top \bigg] \in \mathbb{R}^{(L+K) \times (L+K)}\) is positive definite. This is because (i) \(\pi^*_{u,w}(a|s) > 0\) for all \(a\) from 7 (which has this favorable property that facilitates analysis), and (ii) for any \(y \in \mathbb{R}^{K+L} ~ \text{with} ~ y \neq \mathbf{0}\), \(y^\top B^{\pi^*_{u,w}}(s) y = \sum_{a} \pi^*_{u,w}(a|s) \bigg( y^\top ( Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot|s)} [Q^{\pi^*_{u,w}}(s,a')] ) \bigg)^2 > 0\) as at least one \(a\) should satisfy \(y^\top ( Q^{\pi^*_{u,w}}(s,a) - \mathbb{E}_{a' \sim \pi^*_{u,w}(\cdot|s)} [Q^{\pi^*_{u,w}}(s,a')] ) \neq 0\).

By Proposition 5, we have the Hessian of \(\mathcal{L}(u,w)\) as \(H[\mathcal{L}(u,w)] = \frac{1}{\beta} \sum_{l=1}^p [ \mu_0^\top(I_p - \gamma T^{\pi^*_{u,w}})^{-1}]_l B^{\pi^*_{u,w}}(s_l) = \frac{1}{\beta} \sum_s \rho^{\pi^*_{u,w}}(s) B^{\pi^*_{u,w}}(s)\) where \(p = |\mathcal{S}|\) and \(\rho^{\pi^*_{u,w}}(s) = \sum_{t=0}^\infty \gamma^t \text{Pr}(s_t=s | \pi^*_{u,w}, T, \mu_0)\). Therefore, \(H[\mathcal{L}(u,w)]\) is positive definite under the assumption.

12.2 Existence of \(\lambda\)↩︎

Let \(R = \{ \rho \in \mathbb{R}^{|\mathcal{S}| |\mathcal{A}| }_+ | \sum_{a'} \rho(s',a') = \mu_0(s') + \gamma \sum_{(s,a)} T(s' | s,a) \rho(s,a) ,~\forall s' \}\) be the set of occupancy measures, which is convex and closed. The set \(R\) is also bounded since \(\rho(s,a) \leq \frac{1}{1-\gamma} ~ \forall (s,a)\) for \(\gamma \in [0,1)\). Under Slater condition with strict feasibility, there exists an occupancy measure \(\bar{\rho} \in R\) and positive constants of \(\{ \nu_l \}_{l=1}^L\) such that \[\sum_{(s,a)} c^{(l)}(s,a) \bar{\rho}(s,a) = J_{K+l}(\pi^{\bar{\rho}}) \geq C^{(l)} + \nu_l, ~~ l =1, \cdots, L.\] Although \(\sum_{(s,a)} c^{(l)}(s,a)\rho(s,a)\) is linear in \(\rho\), so that satisfying equality already implies Slater condition, assuming strict feasibility is a mild condition and is standard in prior constrained MDP literature [29].

Define the Lagrangian of the primal problem as: \[\tilde{L}(\rho, u) = \min_{1 \leq k \leq K} \bigg( \sum_{(s,a)} r^{(k)}(s,a) \rho(s,a) \bigg) + \beta \sum_{s} \mathcal{H}_\rho(s) \rho(s) + \sum_{l=1}^L u_l ( \sum_{(s,a)} c^{(l)}(s,a) \rho(s,a) - C^{(l)} )\] with \(\rho \in R\) and \(u \in \mathbb{R}^L_+\). Let the Lagrangian dual be \(D(u) = \max_{\rho \in R} \tilde{L}(\rho, u)\). Let \(P^* = \min_{1 \leq k \leq K} \bigg( \sum_{(s,a)} r^{(k)}(s,a) \rho^*(s,a) \bigg) + \beta \sum_{s} \mathcal{H}_{\rho^*}(s) \rho^*(s)\) denote the optimal value of the primal problem, where \(\rho^* \in R\) is an optimal solution that satisfies the \(L\) constraints.

By strong duality, we have \(P^* = \min_{u \in \mathbb{R}^L_+} D(u) = \tilde{L}(\rho^*, u^*)\) where \(u^* \in \mathbb{R}^L_+\) is an optimal dual variable corresponding to the constraints: \(\sum_{(s,a)} c^{(l)}(s,a) \rho(s,a) \geq C^{(l)}, ~ l =1, \cdots, L\). Also, \[P^* = \max_{\rho \in R} \tilde{L}(\rho, u^*) \geq \tilde{L}(\bar{\rho}, u^*) \geq \min_{1 \leq k \leq K} \bigg( \sum_{(s,a)} r^{(k)}(s,a) \bar{\rho}(s,a) \bigg) + \beta \sum_{s} \mathcal{H}_{\bar{\rho}}(s) \bar{\rho}(s) + \sum_{l=1}^L u_l^* \nu_l\] where the second inequality follows from the definition of \(\{ \nu_l \}_{l=1}^L\). Thus, \[0 \leq u_l^* \leq \frac{P^* - (\min_{1 \leq k \leq K} J_k(\pi^{\bar{\rho}}) + \beta \sum_{s} \mathcal{H}_{\bar{\rho}}(s) \bar{\rho}(s))}{\nu_l} := \bar{u}_l, \forall l.\]

Therefore, it suffices to restrict attention to the compact set \(u \in \Pi_{l=1}^L [0, \bar{u}_l] =: \mathcal{C}\) instead of \(\mathbb{R}^L_+\) when searching for an optimal dual variable corresponding to the constraints \(\sum_{(s,a)} c^{(l)}(s,a) \rho(s,a) \geq C^{(l)}, ~ l =1, \cdots, L\).

Next, the Hessian is explicitly written as follows: \[H[\mathcal{L}(u,w)] = \frac{1}{\beta} \sum_{l=1}^p [ \mu_0^\top (I - \gamma T^{\pi^*_{u,w}})^{-1} ]_{l} B^{\pi^*_{u,w}}(s_l).\]

From the differentiability of \(v^*_{u,w}\) established in Theorem 3.3, \(v^*_{u,w}\) is continuous in \((u,w)\). By their definitions, this continuity extends to \(Q^*_{u,w}\), \(\pi^*_{u,w}\), each entry of \(B^{\pi^*_{u,w}}\) and \(H[\mathcal{L}(u,w)]\). For two real symmetric matrices \(A\) and \(B\), the minimum eigenvalue \(\lambda_{\text{min}}(\cdot)\) satisfies: \[|\lambda_{\text{min}}(A) - \lambda_{\text{min}}(B)| \leq \| A - B \|_2 \leq \| A - B \|_F,\] where \(\| \cdot \|_F\) is the Frobenius norm.

Therefore, \(\lambda_{\text{min}}(H[\mathcal{L}(u,w)])\) is continuous in \((u,w)\) over the compact set \(\mathcal{C} \times \Delta^K\). By the extreme value theorem, there exists \(\lambda = \min_{(u,w) \in \mathcal{C} \times \Delta^K} \lambda_{\text{min}}(H[\mathcal{L}(u,w)]) > 0\) (with \(\lambda \leq \alpha\)).

12.3 Proof of Convergence Analysis↩︎

Let \(\lambda_{\text{min}}( A )\) be the minimum eigenvalue of a real symmetric matrix \(A\). For simplicity, we denote \(\lambda := \lambda_{\text{min}}( H[\mathcal{L}(u,w)] )\). Then \(0 < \lambda \leq \alpha\) [35] and \(\mathcal{L}(u,w)\) is \(\lambda\)-strongly convex.

Theorem 6 Let \((u^*,w^*)\) denote an optimal solution to ?? . For each outer-loop index \(m \geq 1\) in Algorithm 1, let \(Q^*_{u^m, w^m}\) denote the fixed point of 10 with \([u;w] = [u^m;w^m]\), and let \(Q^m\) denote the Q-function after completing the \(m\)-th inner-loop update. For each \(m\), assume \(\| Q^m - Q^*_{u^m, w^m} \|_\infty < \epsilon\) for some \(\epsilon > 0\). Then for \(m \geq 1\), \[\| [u^m;w^m] - [u^*;w^*] \|_{2} \leq (1-\frac{\lambda}{\alpha})^m \| [u^0;w^0] - [u^*;w^*] \|_{2} + \frac{\sqrt{|\mathcal{S}|}}{\lambda} \sqrt{\sum_{i=1}^{K+L} \{ r^{(i)}_{\text{max}} \}^2} \frac{1+\gamma }{(1-\gamma)^2} \epsilon.\]

12.4 Convergence Analysis for Degenerate Case↩︎

We note that the error of \(\mathcal{L}( \frac{1}{m} \sum_{i=1}^m (u^i, w^i) ) - \mathcal{L}(u^*,w^*)\) decreases at rate \(O(\frac{1}{m})\) when \(\{\epsilon_i\}_{i=1}^\infty\) is summable (e.g., \(\epsilon_m = O(\frac{1}{m^{1+\delta}})\) with \(\delta > 0\)).

13 Experimental Details: Tabular Settings↩︎

When generating structured MOMDPs randomly, we first verify whether the generated instances are feasible. To do this, We first consider the following unregularized convex optimization: \[\label{eq:optimal95obj} \max_{\rho \geq 0} \min_{1 \leq k \leq K} \bigg( \sum_{(s,a)} r^{(k)}(s,a) \rho(s,a) \bigg)\tag{32}\] \[\label{eq:optimal95flow} \sum_{a'} \rho(s',a') = \mu_0(s') + \gamma \sum_{(s,a)} T(s' | s,a) \rho(s,a) ,~\forall s'\tag{33}\] \[\label{eq:optimal95const95c} \sum_{(s,a)} c^{(l)}(s,a) \rho(s,a) \geq C^{(l)}, ~~ l =1, \cdots, L\tag{34}\] which is equivalently expressed as the following LP by using additional scalar variable \(\tilde{c} \in \mathbb{R}\): \[\label{eq:eqiv95optimal95obj} \max_{\rho \geq 0, \tilde{c}} \tilde{c}\tag{35}\] \[\label{eq:eqiv95optimal95flow} \sum_{a'} \rho(s',a') = \mu_0(s') + \gamma \sum_{(s,a)} T(s' | s,a) \rho(s,a) ,~\forall s'\tag{36}\] \[\label{eq:eqiv95optimal95const95r} \sum_{(s,a)} r^{(k)}(s,a) \rho(s,a) \geq \tilde{c}, ~~ k =1, \cdots, K,\tag{37}\] \[\label{eq:eqiv95optimal95const95c} \sum_{(s,a)} c^{(l)}(s,a) \rho(s,a) \geq C^{(l)}, ~~ l =1, \cdots, L.\tag{38}\]

We want to generate \(\mu_0, T, r\), and \(c\) in structured MOMDPs to satisfy feasibility and Slater condition by solving the following LP using the pywraplp function from the OR-Tools library:

\[\label{eq:feasible95obj} \max_{\rho \geq \epsilon_{\text{low}}} ~ 0\tag{39}\] \[\label{eq:feasible95flow} \sum_{a'} \rho(s',a') = \mu_0(s') + \gamma \sum_{(s,a)} T(s' | s,a) \rho(s,a) ,~\forall s'\tag{40}\] \[\label{eq:feasible95const95r} \sum_{(s,a)} r^{(k)}(s,a) \rho(s,a) \geq \tilde{c} + \epsilon_{\text{low}}, ~~ k =1, \cdots, K,\tag{41}\] \[\label{eq:feasible95const95c} \sum_{(s,a)} c^{(l)}(s,a) \rho(s,a) \geq C^{(l)} + \epsilon_{\text{low}}, ~~ l =1, \cdots, L\tag{42}\]

where \(\epsilon_{\text{low}}\) is used to guarantee the strict feasibility for Slater condition, and we set \(\epsilon_{\text{low}}=10^{-4}\). If the LP solver does not find a feasible solution, we regenerate the constrained MOMDP until a feasible instance is found. Once any feasible solution is found, we solve the LP of 35 , 36 , 37 , and 38 by using LP solver to acquire the optimal max-min value \(\tilde{c}^*\).

Each algorithm is updated iteratively until the maximum change in the \(Q\)-function between successive iterations falls below \(10^{-4}\). We use the following settings: \(\gamma=0.8\), \(l_{\text{w}}=0.001\), and \(\text{ITER}=3000\). For each algorithm, we independently disable the learning of \(u\), of \(w\), and of both \((u,w)\), initializing \(u\) to the zero vector and \(w=[1/K, \cdots, 1/K] \in \Delta^K\). Each algorithm is evaluated using three random seeds for each constrained MOMDP setting, resulting in six runs when averaged across the MOMDP classes. All experiments are conducted on an Intel Core i9-10900X CPU @ 3.70GHz.

14 Experimental Details: Applications↩︎

14.1 Implementation of Our Algorithm for Applications↩︎

We now leverage the usage of neural network for our algorithm. If we differentiate the both side of \(v^{*}_{u, w}(s) = [\mathcal{T}_{u,w} v^{*}_{u, w}](s)\) w.r.t. \(u\) and \(w\) for all \(s\), then we have the following formula: \[\label{eq:grad95w95v95star} \forall s, ~ \nabla_w v^{*}_{u, w}(s) = \sum_a \pi^*_{u,w}(a | s) \bigg( r(s,a) + \gamma \sum_{s'}T(s'|s,a) \nabla_w v^{*}_{u, w}(s') \bigg).\tag{43}\]

\[\label{eq:grad95u95v95star} \forall s, ~ \nabla_u v^{*}_{u, w}(s) = \sum_a \pi^*_{u,w}(a | s) \bigg( c(s,a) + \gamma \sum_{s'}T(s'|s,a) \nabla_u v^{*}_{u, w}(s') \bigg).\tag{44}\]

Here, \(\pi^*_{u,w}(a | s)\) is defined as in 7 . To ensure stable gradient estimation in continuous state spaces, we parameterize a gradient network to estimate \(\nabla_u v^{*}_{u, w}(s)\) and \(\nabla_w v^{*}_{u, w}(s)\). For continuous action environments, we also employ an actor network \(\pi_\theta\) and implement Algorithm 2. To further stabilize the estimation of the gradient, we add an additional linear layer after the penultimate layer of the actor network \(\pi_\theta\), and use its \((L+K)\)-dimensional output as the gradient network \(g_\theta(s)\). We use the notation \(g_\theta\) to indicate that the actor network and the gradient network share parameters and jointly update their lower-layer weights.

14.2 Unconstrained Max-min MORL Algorithm↩︎

14.3 Building Thermal Control↩︎

We use a building thermal control environment from SustainGym [54], which simulates a commercial building with \(N_{\text{zone}}=23\) thermal zones. We group these zones into three sections (bottom, middle, and top) and define a comfort objective for each section. Accordingly, in our experiments we set \(K=3\) and \(L=1\), i.e., three reward objectives and a single cost constraint.

At each timestep \(t\), the agent observes a normalized continuous state \(s_t \in \mathbb{R}^{29}\). The underlying (unnormalized) state is constructed by concatenating: (i) the temperatures of all 23 zones, (ii) outdoor temperature, (iii) global horizontal irradiance (GHI), (iv) ground temperature, (v) occupancy power, (vi) carbon intensity, and (vii) electricity price. The environment clips each component to predefined bounds and linearly rescales it to \([0,1]\) before returning \(s_t\).

The action is a continuous vector \(a_t \in \mathbb{R}^{23}\) with each component bounded in \([-1,1]\), representing zone-level heating/cooling power commands (cooling is negative and heating is positive). The reward is a \(K\)-dimensional vector \(\mathbf{r}_t \in \mathbb{R}^{3}\) that measures comfort for three zone groups. The 23 zones are partitioned into bottom (zones 1–9), middle (zones 10–16), and top (zones 17–23), and each reward entry increases when the temperatures in the corresponding group stay closer to the target temperatures. We define a single scalar energy cost signal \(c_t\) from the aggregated magnitude of the zone-level power commands, specifically \(c_t = \|a_t\|_1 / s\), where we tune the scaling factor \(s\) (set to \(s=5\)) to match the reward scale, and use \(c_t\) as the cost for the energy constraint (\(L=1\)).

Each episode consists of 480 timesteps. The total training budget is \(8\times 10^{5}\) timesteps, with evaluations conducted once every two training episodes. During each evaluation, one episode is run and the cumulative discounted sums of the \((L+K)\)-dimensional vector are computed. These experiments were conducted using an NVIDIA RTX 4000 SFF Ada GPU (20GB) across five random seeds.

14.4 Effect of \(\beta\) and Wall-Clock Time Comparison↩︎

The entropy term is used not only to encourage exploration but also to resolve the indeterminacy issue in solving the dual of the max-min optimization problem [12]. The impact of \(\beta\) is seen in Table 8: The max-min performance in the Building environment is higher for moderate values of \(\beta\) (e.g., 0.1, 0.05) than for very small values (e.g., \(\beta=0.01\)).

Table 9 reports the wall-clock training time per episode for each algorithm. Incorporating constraints increases the training time from MA-SAC to MA-SAC-L. Although our algorithm requires more computation than MA-SAC-L, it achieves a better balance between max-min fairness and constraint satisfaction, as demonstrated in Table 3. Furthermore, our algorithm requires significantly less computation than Max-min GS.

14.5 Multi-Objective Locomotion Control↩︎

14.6 Evaluation with Different \(C_{th}\)↩︎

We set the threshold value \(C_{th} = 40\) in the MoAnt environment, which is more stringent than the value \(C_{th} = 50\) used in Table 5. Our method consistently balances max-min fairness with constraint satisfaction, meeting the constraint while achieving superior max-min performance compared to MA-SAC-L.

14.7 Greenhouse-gas-emission-aware Traffic Management↩︎

Here, \(|\mathcal{A}|=4\). MA-PGO is implemented using the PPO algorithm [64] to maximize the average reward over \(K=16\) objectives. Based on MA-PGO, MA-CPGO applies clipping only to the unconstrained reward part (i.e., \(w \cdot r\)), while leaving the constrained reward part unclipped (i.e., \(u \cdot c\)) and applies the Lagrangian update, following [25] to improve constraint satisfaction. We found that applying clipping to \(u \cdot c + w \cdot r\) does not guarantee constraint satisfaction. Each method is run for 100k timesteps per seed using five random seeds.

15 Limitation, Future Work, and Discussion↩︎

In this section, we discuss several limitations of our work and related future research avenues, although our method offers a promising direction for developing constrained MORL algorithms.

First, there is a lack of well-established benchmarks for MORL compared to standard RL settings [3], and even fewer environments are specifically designed for constrained MORL. Additionally, most existing MORL environments have low-dimensional reward spaces (typically fewer than four dimensions) [42], which limits the ability to evaluate our algorithm in high-dimensional settings. Developing practical benchmarks for both MORL and constrained MORL is therefore a critical research direction for the community.

Second, while it is common in the constrained MDP literature to assume that feasibility is ensured by appropriately chosen thresholds [26], [31], determining such thresholds, that is, setting the constraint set \(\{ C^{(l)} \}_{l=1}^L\), is non-trivial in practice outside of simple or tabular domains. Unlike trial-and-error reward design, constraint threshold design is often infeasible or unsafe due to the potential risks and costs involved. Leveraging external sources of information, such as human demonstrations or natural language descriptions, offers a promising path for setting constraint thresholds in constrained RL and MORL. Another possible approach is to infer the constraint values from expert demonstrations, commonly referred to as inverse constrained RL [65], [66].

Third, while our setting clearly distinguishes rewards from costs, this distinction may be ambiguous in other domains. Determining which objectives should be treated as constraints versus unconstrained rewards can be challenging. As with constraint threshold design, incorporating external guidance could help better structure constrained MORL problems.

Fourth, several constrained RL studies have explored more conservative formulations than those based on expected cumulative cost, for example, using outage probability or quantile-based constraints to manage rare but critical failures in domains such as finance or insurance [56], [67]. While our current framework and analysis rely on expected cumulative cost, extending it to support such conservative constraint formulations presents a valuable direction for future work.

Lastly, although we assume the convergence of the (action) value function for each weight pair \((u,w)\), it is well known that the combination of function approximation, bootstrapped updates, and off-policy learning can lead to instability and even divergence during training [61], [68]. A theoretical investigation into this so-called deadly triad, along with additional convergence guarantees, would further improve the robustness of our algorithm and broaden its applicability to other domains.

Scalarization-based methods are highly valuable, especially because of their interpretability and flexibility in expressing designer preferences. In particular, when incorporating constraints, these methods also make it straightforward to assess constraint satisfaction through the corresponding dual variables. However, linear scalarization cannot recover nonlinear or concave regions of the Pareto frontier, potentially missing desirable trade-off solutions [1], [3]. While mixtures of convex scalarization functions can help approximate concave regions, this often requires careful tuning and may increase computational effort. Addressing the limitations of scalarization-based approaches is indeed valuable.

Finally, we note that the purpose of our theoretical analysis is not to imply that all assumptions will be checked analytically in practice, but rather to provide predictable behavior and guidance for practical usage of our algorithm.

References↩︎