Model-Based Reinforcement Learning for Infinite-Horizon Discounted Constrained Markov Decision Processes

1 Introduction

Markov Decision Processes (MDPs) are a powerful approach to model stochastic systems where stationary control policies are appropriate. Many cyber-physical (i.e. physical systems controlled algorithmically) systems bear intrinsic limitations on the nature of control that may be applied. Hence, Constrained-MDP (CMDP) are an appropriate framework for the modeling and analysis of such systems [Altman, 1999].

In this paper, we aim to develop simple algorithms to learn near-optimal policies for a CMDP without knowing the system parameters. Although, a regular model-based RL algorithm attempts to collect as few samples as possible to quickly solve for the optimal policy, minimizing the number of samples taken is even more essential in the CMDP setting. This requirement is due to the existence of constraints in the CMDP setting, and it might be important to violate them as few times as possible while maximizing the objective of the system. Therefore, the behavior of a system with respect to (w.r.t) both objective maximization and safety violation over time is a crucial performance metric for a proposed RL algorithm for CMDPs.

Main Contributions: Our goal is to upper bound the number of samples required to learn a near-optimal policy while nearly satisfying the constraints with high probability (w.h.p.) in the context of the discounted infinite-horizon setting. Our contributions are mainly threefold:

• (i) We design and analyze two model-based RL algorithms for CMDPs. One of them pursues a generative model based approach that obtains samples initially and creates a model. The other one is based on an online approach in which the model is updated over time-steps. With both algorithms, the estimated model might lead to infeasible situation. Thus, we utilize the idea of a confidence-ball around the estimated model such that the true model would belong to that ball w.h.p. This ensures that a solution may be found w.h.p. under the assumption that the real model has a solution.
• (ii) Both algorithms follow a two-stage pattern of model construction and a CMDP solution. The algorithms use linear programming (LP) to solve the CMDP problem with additional linear constraints to incorporate the confidence-ball.

(iii) We characterize PAC-type sample complexity bounds for both algorithms, accounting for both objective maximization and constraint satisfaction. Intuitively, the model constructed by these algorithms must be more accurate than models created by unconstrained counterparts, which conjecture our main results are consistent with. Furthermore, a comparison of our main findings with lower bounds on sample complexity of MDPs [Azar et al., 2013; Dann and Brunskill, 2015] shows an increase in our results by a logarithmic factor in the number of constraints and the size of the state space. However, there is no earlier work on lower bound of sample complexity of learning CMDPs to our best knowledge.

As mentioned above, cyber-physical systems might have a large number of constraints. However, our results indicate that the number of constraints should not be a major concern in implementation, since our bounds scale logarithmically with number of constraints. Hence, the results suggest that the constrained RL approach is likely applicable in a straightforward manner to cyber-physical systems.

RelatedWork: There are many articles studying the problem of controlling CMDPs with an algorithmic approach and control-theoretic view [Altman, 1999; Altman, 2002; Borkar, 2005; Borkar and Jain, 2014; Singh and Kumar, 2018; Singh et al., 2014]. The results take the form of proving asymptotic convergence of their proposed methods under the assumption of the known model. There are also extensions of this approach to the context of an unknown model, where the focus is still on asymptotic behavior [Bhatnagar and Lakshmanan, 2012; Chow et al., 2018; Tessler et al., 2018; Paternain et al., 2019]. These studies use Lagrangian method to show zero duality gap asymptotically. Further, [Liu et al., 2019] also develops an algorithm based on the Lagrangian method, but with small eventual duality gap. Finally, empirical studies based on the Lagrangian method have also been presented [Liang et al., 2018].

There are also studies on the constrained bandit case. Although bandits are not MDPs per se, they are strongly related to them. Articles such as [Badanidiyuru et al., 2013; Wu et al., 2015; Amani et al., 2019] consider such constraints, either in a knapsack sense, or on the type of controls that may be applied in a linear bandit context.

More related to our work theme are parallel studies on CMDPs. For example, [Zheng and Ratliff, 2020] and [Wachi and Sui, 2020] provide results with the assumption of unknown reward functions, with either a known or deterministic transition kernel. There are other works [Satija et al., 2020] focusing on proving asymptotic convergence without providing a bound on learning rate. Finally, closest related work to this article is [Efroni et al., 2020] which explores algorithms similar to ours in finite-horizon setting, but concentrating on characterizing objective and constrained regret bounds. Now, regret and sample complexity bounds are not directly translatable [Dann et al., 2017], and converting their regret bounds to our setting gives relatively weak sample complexity bounds. Specifically, our main results with logarithmic increase in sample complexity with the number of constraints differentiates our work.

2 Notation and ProblemFormulation

Notation and Setup: Our focus is on an infinite-horizon CMDP defined by a tuple M = 〈S,A, P, r, c, ¯ C, s0, γ〉. S and A represent the sets of states and actions respectively. Additionally, P(s′ |s, a) is used to indicate the probability of reaching state s′ by taking action a while being at state s. We define r(s, a) as the reward for each state-action pair (s, a). We assume that there are N constraints. We use c to denote the cost matrix, where c(i, s, a) is the immediate cost incurred by the ith constraint in (s, a) where i ∈ {1, . . . , N}. Further, the value of the constraints (i.e. the bound that must be satisfied) are determined by the vector ¯ C. Also, initial state is specified by s0. Finally, we use γ for discount factor. In this study, the discount factor is unique for both objective function and constraint functions where they shall be defined later.

Assumption 1. State and action sets S and A are assumed to be finite with cardinalities |S| and |A|. In addition, the immediate cost and immediate reward r(s, a) are assumed to be taken from the interval [0, 1]. Number of constraints is also assumed to be N which for each i ∈ {1, . . . ,N}, ¯ Ci ∈ [0, ¯ Cmax].

Now, we define a stationary policy π : S×A → [0, 1]|A| as a mapping from state-action space S ×A to set of probability vectors defined over action space in order to choose an action at any time-step t. Henceforth, π(s, a) represents the probability of choosing the action a when the system is at state s. Also, a ∼ π(s, ·) means that action a is chosen according to stationary policy π while being at state s.

Fixing a policy π transforms the underlying MDP to a Markov chain. The transition kernel of this Markov chain is Pπ, which can be viewed as an operator. The operator Pπf(s) = E[f(st+1)|st = s] = ∑s′∈S Pπ(s′ |s)f(s′ ) takes any function f : S → R and returns the expected value of f in the next time-step. For convenience, we define the multistep version P tπ f(s) = PπPπ . . . Pπf, which is repeated t times. Further, we define P 0 π as the identity operator.

For the objective and constraint functions, we consider discounted infinite-horizon criteria with identical discount factor γ. We define the value function of state s under policy π as

where expectation E[·] is taken w.r.t transition kernel P. Next, the local variance of the value function at time step t under

Analogous to the definition of the value function (1), the ith constraint function under policy π is defined as

Again, the local variance of ith constraint function under policy π, i.e. σ2Cπ i is defined similar to local variance of value function (2).

Assumption 2. We assume that the CMDP problem of (4) is feasible with optimal policy π∗ and optimal solution V ∗ (s0) = V π∗ (s0).

Note that we only consider learning feasible CMDPs by this assumption.

Constrained-RL Problem: The Constrained RL problem formulation is identical to the CMDP optimization problem of (4) with one difference. Here, we are not aware of the values of the transition kernel P.1 We desire to provide modelbased algorithms and determine the sample complexity results in a PAC sense, which is defined as follows:

Definition 1. For an algorithm A, sample complexity is the number of samples that A requires to achieve

for a given ϵ and δ.

1 We only assume that transition kernel is unknown and the extension to unknown reward and cost matrices is straightforward, and does not require additional methodology.

Note that with this definition, we include both objective maximization and constraint violations as opposed to the traditional definition that only considers the objective [Strehl and Littman, 2008].

3 Sample Complexity Result of Generative Model Based Learning

Generative model based learning is a well known approach to learn an optimal policy for an MDP. However, naive application of this approach to CMDPs may not end with a feasible solution. Hence, we explore the generative model based approach for CMDPs, and propose a generative model based CMDP learning algorithm called Generative ModelConstrained RL (GM-CRL). According to GM-CRL, each state-action pair is sampled n number of times uniformly across all state-action pairs, the number of times each transition occurs n(s′ , s, a) for each next state s′ is counted, and an empirical model of transition kernel denoted by ̂ P (s′ |s, a) = n(s′ ,s,a) n ∀(s′ , s, a) is constructed.

Unlike MDP problem formulation, there is no guarantee such that CMDP problem formulation w.r.t. ̂ P is feasible. In order to resolve the feasibility concern, we expand the space of transition kernels to include the true transition kernel P, noting that the CMDP problem w.r.t. P is feasible from Assumption 2. The algorithmic layout of this approach is as follows. GM-CRL creates a class of CMDPs using the empirical model. This class is denoted by MδP and contains CMDPs with identical reward, cost matrices, ¯ C, initial state s0 and discount factor of the true CMDP, but with transition kernels close to true model. This class of CMDPs is defined

where δP is defined in Algorithm 1. Note that for any M ′ ∈ M, objective function V ′ π (s0) and cost functions C ′ π i (s0) are computed w.r.t. the corresponding transition kernel P ′ according to equations (1) and (3) respectively.

At the end, GM-CRL maximizes the objective function among all possible transition kernels, while satisfying constraints (if feasible). More specifically, it solves the optimistic planning problem below

To solve the problem of (7), GM-CRL uses Extended Linear Programming, or ELP. This method takes MδP as input and gives ˜π for the optimal solution. The description of ELP is provided in supplementary materials. Algorithm 1 describes GM-CRL.

Algorithm 1 GM-CRL

• 1: Input: accuracy ϵ and failure tolerance δ.
• 2: Set δP = δ 5(N+2)|S|3 |A| .
• 3: Set n(s′ , s, a) = 0 ∀(s, a, s′ ).
• 4: for each (s, a) ∈ S ×A do
• 5: Sample (s, a), n = 1152(log 2)2 γ 2 ϵ2 (1-γ)3 |S|2 |A| log 4 δP and update n(s′ , s, a).
• 6: ̂ P(s′ |s, a) = n(s′ ,s,a) n ∀s′ .
• 7: end for
• 8: ConstructMδP according to (5).
• 9: Output ˜π = ELP(MδP ).

3.1 PACAnalysis of GM-CRL

Here, we present the sample complexity result of GM-CRL.

Theorem 1. Consider any infinite-horizon CMDP M = 〈S,A, P, r, c, ¯ C, s0, γ〉 satisfying assumptions 1 and 2, and CMDP problem formulation of (4). Then, for any ϵ ∈ (0, 0.22γ √|S|(1-γ) ) and δ ∈ (0, 1), algorithm 1 creates a model

C, s0, γ〉 and outputs policy ˜π

CMDP ˜M = 〈S,A, ˜ P, r, c, ¯ such that

with at least total sampling budget of

The proof of Theorem 1 is different from the traditional analysis framework of unconstrained RL [Azar et al., 2013] in the following manner. First, consider the role played by optimism in model construction. The notion of optimism is not required for learning unconstrained MDPs with generative models, because any estimated model is always feasible [Puterman, 2014]. However, there is no such guarantee for a general CMDP problem formulation [Altman, 1999]. Specifically, simply substituting the true kernel P by the estimated one ̂ P is not appropriate, since there is no assurance of feasibility of that problem. Hence, GM-CRL converts the CMDP problem under the estimated transition kernel to an optimistic planning problem (7) and an ELP-based solution.

Second, the core of the analysis of every unconstrained MDP is based on being able to characterize the optimal policy via the Bellman operator. This technique enables one to obtain a sample complexity that scales with the size of the state space as O(|S|). However, we cannot use this approach to characterize the optimal policy in a CMDP [Altman, 1999]. We require a uniform PAC result over set of all policies and set of value and constraint functions, which in turn leads to quadratic sample complexity in the size of state space; i.e., a scaling of O(|S|2 ).

Corollary 1. In case of N = 0, the problem would become regular unconstrained MDP. And, the sample complexity result with N = 0 would also hold for unconstrained case.

Now, we present some of the lemmas that are essential to prove Theorem 1. Then we sketch the proof of this theorem. The detailed proofs are provided in supplementary materials.

First, we show that true CMDP lies inside the MδP with high probability, w.h.p. Hence, the problem (7) is feasible w.h.p., since the original CMDP problem is assumed to be feasible according to Assumption 2.

Lemma 1.

Proof Sketch: Fix a state-action pair (s, a) and next state s′ . Then, according to combination of Hoefding's inequality [Hoeffding, 1994] and empirical Bernstein's inequality [Maurer and Pontil, 2009], we obtain that each P(s′ |s, a) is inside the confidence set defined by (6) with probability at least 1- δP . Applying the union bound yields the result. ✷

Now, we present the core lemma required for proving Theorem 1 and its proof sketch. Using this lemma, we bound the mismatch in objective and constraint functions when we have n number of samples from each (s, a). This bound applies uniformly over the set of policies and set of value and constraint functions. The result also enables us to bound the objective and constraint functions individually. Then we apply the union bound on all objective and constraint functions. This process is the reason why the number of constraints appear logarithmically in the sample complexity result.

Lemma 2. Let δP ∈ (0, 1). Then, if n ≥ 11819 |S|2 log 4/δP (1-γ)2 , under any policy π

w.p. at least 1-5|S|3 |A|δP .

Proof Sketch: We first show that | ˜ P(s′ |s, a) - P(s′ |s, a)| ≤ O(√ P(s′ |s,a)(1-P(s′ |s,a)) n ) for each s′ , s, a. Then, we show that (Pπ - ˜ Pπ)V π (s) ≤ O(√ |S| n σV π(s)). Applying this bound to | ˜ V π (s0) - V π (s0)| and from the fact that σV π(s) is close to ˜σV π(s) by O( √|S| (1-γ)n1/4 ), we obtain the result. This procedure is also applicable to each constraint function i. ✷

Proof Sketch of Theorem 1: From Lemma 1, we know that the optimistic planning problem (7) is feasible w.h.p. Hence, we can obtain an optimistic policy ˜π. The rest of this proof consists of two major parts.

First, we prove ϵ-optimality of objective function w.h.p. Considering policy π∗ we obtain |V π∗ (s0) - ˜ V π∗ (s0)| ≤ O(√ |S| (1-γ)3 n ) w.h.p. by means of Lemma 2. Similarly,

|V ˜π (s0) - ˜ V ˜π (s0)| ≤ O(√ |S| (1-γ)3 n ) w.h.p. Next, we use the fact that ˜ V π∗ (s0) ≤ ˜ V ˜π (s0) and obtain

Next, we show that each constraint is violated at most by ϵ w.h.p. Here, we use the second part of Lemma 2 to bound constraint violation. Thus, for each i ∈ {1, . . . ,N} we have |C ˜π i (s0) - ˜ C ˜π i (s0)| ≤ O(√ |S| (1-γ)3 n ) w.h.p. Also, we know that ˜ C ˜π i (s0) ≤ ¯ Ci, since ˜π is solution of the ELP. Hence, we obtain

w.h.p. Finally, we obtain the end result by applying the union bound, and obtaining n by solving ϵ = O(√ |S| (1-γ)3 n ). ✷

4 Sample Complexity Result of Online Learning

The GM-CRL approach operates in a way that every stateaction pair in the system is sampled a certain number of times before a policy is computed. However, there are applications that are not capable of utilizing this approach, since it may not be possible to reach those states without the employment of some policy, or they might be unsafe, and so should not be sampled often. Hence, we have to find an approach that can collect samples from the environment by means of an online algorithm.

Upper Confidence Constrained-RL, or UC-CRL described in Algorithm 2, is an online method proceeding over timesteps. At each time-step t, UC-CRL constructs an empirical model ̂ P using state-action visitation frequencies, i.e., ̂ P(s′ |s, a) = n(s′ ,s,a) n(s,a) , where n(s′ , s, a) and n(s, a) are visitation frequencies. Then, we use ̂ P to create a confidence interval around each element ̂ P(s′ |s, a) using same concentration inequalities of GM-CRL defined by (6). Next, UCCRL constructs set of infinite-horizon CMDPs Mt which any CMDP M ′ ∈ Mt has identical discount factor and reward and cost matrices to the true CMDP M, but different transition kernels from the concentration inequalities. Mt is identical toMδP except for the use of n(s, a) instead of n. Thus the class of CMPDs is defined as below at each time-step t :

where δ1 is defined in Algorithm 2.

Subsequently, UC-CRL uses ELP to solve the optimistic CMDP problem below and get the optimistic policy ˜πt :

This problem is identical to the problem of (7), except for substituting MδP with Mt. Here, for any M ′ ∈ Mt, V ′π (s0) and C ′π i (s0) are computed according to (1) and (3) w.r.t. underlying transition kernel P ′ , respectively.

Algorithm 2 UC-CRL

• 1: Input: accuracy ϵ and failure tolerance δ.
• 2: Set m according to (9) and (10).
• 3: Set t = 1,wmin = ϵ(1-γ) 4|S| ,Umax = |S|2 |A|m, δ1 = δ 4(N+1)|S|Umax .
• 4: Set n(s, a) = n(s′ , s, a) = 0 ∀s, s′ ∈ S, a ∈ A.
• 5: while there is (s, a) with n(s, a) < |S|m 1-γ do
• 6: ̂ P(s′ |s, a) = n(s′ ,s,a) n(s,a) ∀(s, a) with n(s, a) > 0 and s′ ∈ S.
• 7: ConstructMt according to (8).
• 8: ˜πt = ELP(Mt).
• 9: at ∼ ˜πt(st), st+1 ∼ P(·|st, at)
• 10: if n(st, at) < |S|m 1-γ then
• 11: n(st, at) + +, n(st+1, st, at) + +.
• 12: end if
• 13: t ++
• 14: end while

UC-CRL is inspired by the infinite-horizon algorithm UCRL-γ [Lattimore and Hutter, 2014] and its finite-horizon equivalent UCFH [Dann and Brunskill, 2015] with differences. Similar to UCRL-γ, Algorithm 2 uses a combination of the empirical Bernstein's and Hoeffding's inequalities. These concentration inequalities allow us to ensure linearity of constraints (i.e., we can indeed use an extended linear program to solve for the constrained optimistic policy). However, the constraints of UCFH contain non-linear expressions preventing us from employing ELP. Furthermore, unlike UCRL-γ and UCFH, Algorithm 2 updates the model at each time-step rather than at the beginning of long phases. This procedure allows for faster model construction. Finally, since we are solving a CMDP, this algorithm utilizes ELP instead of Extended Value Iteration which is used by UCRL-γ.

4.1 PACAnalysis of UC-CRL

We now present the PAC bound of Algorithm 2.

Theorem 2. Consider CMDP M = 〈S,A, P, r, c, ¯ C, s0, γ〉 satisfying assumptions 1 and 2. For any 0 < ϵ, δ < 1, under UC-CRL we have:

time-steps.

We follow an approach motivated by [Lattimore and Hutter, 2014] and its finite-horizon version [Dann and Brunskill, 2015] to prove Theorem 2. However, there are several differences in our technique, and we need to accommodate the frequent model update in our proof. We will show that, unlike existing approaches, we can update the model at each timestep, without increasing the sample complexity. Thus, we are able to obtain PAC bounds that match the unconstrained case, and only increase by logarithmic factor with the number of constraints.

There are also recent works on characterizing the regret of constrained-RL in a finite-horizon setting [Efroni et al., 2020] with an algorithm similar to Algorithm 2. An important emerging question is whether one can immediately convert these regret results into sample complexity bounds? A naive translation of the regret bounds of [Efroni et al., 2020] would give us a PAC result ˜O( |S|2 |A|H4 ϵ2 ). For comparing finite-horizon setting with infinite-horizon one, we can replace H with 1 1-γ to obtain a PAC result for the equivalent infinite-horizon algorithm. Considering this, the approach followed by [Efroni et al., 2020] gives a PAC bound which is looser than our result by a factor of 1 (1-γ)2 . Therefore, this alternative option does not lead to the strong bounds that we are able to obtain, and matches existing PAC results of the unconstrained case.

Now, we present the notions of knownness and importance for state-action pairs and base our proof on these notions. Then we present the key lemmas needed for proving Theorem 2. Finally, we provide a proof sketch for Theorem 2. The detailed analysis is provided in supplementary materials.

Let the weight of (s, a)-pair under any policy π be its discounted expected frequency

Using this general definition, we define the weight of (s, a) under policy ˜πt aswt(s, a) = w ˜πt (s, a|st).

Then, the importance ιt of (s, a) at time-step t is defined as its relative weight compared to wmin := ϵ(1-γ) 4|S| on a log-scale

Note that ιt(s, a) ∈ {0, 1, 2, 4, 8, 16, . . . } is an integer indicating the influence of the state-action pair on the value function of ˜πt. Similarly, we define knownness as

which indicates how often (s, a) has been observed relative to its importance. Value of m is defined in Algorithm 2. Now, we can categorize (s, a)-pairs into subsets

where Xt = {(s, a) : ιt(s, a) > 0} is the active set and ¯Xt is the set of (s, a)-pairs that are very unlikely under policy ˜πt. We will show that if the criteria |Xt,κ,ι| ≤ κ is met, then the model of UC-CRL would achieve near-optimal policies where these policies would violate constraints at most by ϵ w.h.p. This condition specifies that important state-action pairs under policy ˜πt are visited a sufficiently large number of times. Thus, the model of UC-CRL will be accurate enough to obtain PAC bounds.

Now, we first show that the true model lies inMt for every time-step t w.h.p.

Lemma 3. M ∈ Mt for all time-steps t with probability at least 1- δ 2(N+1) .

Proof Sketch: Let consider a fixed (s, a), next state s′ and a time-step t. Then, P(s′ |s, a) belongs to the confidence set constructed by the combined Bernstein's and Hoeffding's inequalities. By taking the union bound over maximum number of model updates, Umax, and next states we obtain the result. ✷

Next, we bound the number of time-steps in which the condition |Xt,κ,ι| ≤ κ is violated w.h.p.

Lemma 4. Suppose E is the number of time-steps t for which there are κ and ι with |Xt,κ,ι| > κ, i.e. E = ∑∞t=1 I{∃(κ, ι) : |Xt,κ,ι| > κ} and let

where Emax = log2 1 wmin(1-γ) log2 |S|. Then, P(E ≤ 6|S||A|mEmax) ≥ 1- δ 2(N+1) .

Proof sketch: This lemma is proven in two stages. First, we bound the total number of times a fixed (s, a) could be observed in a particular Xt,κ,ι over all time-steps. Then, we provide a high probability bound on the number of time-steps that |Xt,κ,ι| > κ for a fixed (κ, ι). Finally, we get the result using of martingale concentration and union bound. ✷

Finally, the next lemma bounds the mismatch between objective and constraint functions of the optimistic model and true model. This lemma functions similarly to Lemma 2 for GM-CRL. It provides a uniform PAC result over value and constraint functions. Hence, it enables us to have individual PAC results for any objective and constraint functions. As discussed in GM-CRL section, this process is responsible for a logN increase in the PAC result.

Lemma 5. Assume M ∈ Mt. If |Xt,κ,ι| ≤ κ for all (κ, ι) and 0 < ϵ ≤ 1 and

Proof Sketch: First, we show | ˜ P (s′ |s, a) - P (s′ |s, a)| ≤ O(√ P (s′ |s,a)(1-P (s′ |s,a)) n ) for each s′ , s, a. Then we prove that at each time-step t, (Pπ- ˜ Pπ)V π (s) ≤ O(√ |S| n σV π(s)). Next we partition the state-action based on knownness, i.e., whether they belong to Xt or not. By using all bounds and sequence of CMDPs, we obtain a bound on | ˜ V π (s0)-V π (s0)|. Eventually, we use the definition of weights to get the final result. This procedure is also applicable to each constraint ✷

function i.

Proof Sketch of Theorem 2: We first use Lemma 3 and show that true CMDP is admissible ,i.e. M ∈ Mt for every time-step, w.p. at least 1- δ 2(N+1) . Hence, the optimistic planning problem becomes feasible and an optimistic policy ˜πt exists w.h.p. Further, we use Lemma 4 to bound the number of time-steps where |Xt,κ,ι| > κ w.h.p. Thus, for other time-steps where |Xt,κ,ι| ≤ κ, we apply Lemma 5 we show that objective function is ϵ-optimal and all constraint functions are violated by ϵ. Eventually, we obtain the result by means of union bound. ✷

5 Conclusion

In this paper, we presented the notion of sample complexity in objective maximization and constraint satisfaction in order to understand the performance of RL algorithms for safety-constrained applications. We proposed and analyzed two types of algorithms-GM-CRL and UC-CRL. The main finding of a logarithmic factor increase in sample complexity compared to unconstrained regime indicates the value of the algorithms in applying them to real systems.

Acknowledgments

This work was funded in part by grants NSF-Intel CNS 1719384, NSF CPS-2038963, NSF CNS 1955696, ARO W911NF-19-1-0367, NSF CRII- CPS-1850206, NSF CAREER-EPCN- 2045783.