{'fraction_non_alphanumeric': 0.03827934612962405, 'fraction_numerical': 0.01221215191259823, 'mean_word_length': 4.576729271644526, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chance Constrained Policy Optimization for Process Control and Optimization

Chemical process optimization and control are affected by 1) plant-model mismatch, 2) process disturbances, and 3) constraints for safe operation. Reinforcement learning by policy optimization would be a natural way to solve this due to its ability to address stochasticity, plant-model mismatch, and directly account for the effect of future uncertainty and its feedback in a proper closed-loop manner; all without the need of an inner optimization loop. One of the main reasons why reinforcement learning has not been considered for industrial processes (or almost any engineering application) is that it lacks a framework to deal with safety critical constraints. Present algorithms for policy optimization use difficult-to-tune penalty parameters, fail to reliably satisfy state constraints or present guarantees only in expectation. We propose a chance constrained policy optimization (CCPO) algorithm which guarantees the satisfaction of joint chance constraints with a high probability - which is crucial for safety critical tasks. This is achieved by the introduction of constraint tightening (backoffs), which are computed simultaneously with the feedback policy. Backoffs are adjusted with Bayesian optimization using the empirical cumulative distribution function of the probabilistic constraints, and are therefore self-tuned. This results in a general methodology that can be imbued into present policy optimization algorithms to enable them to satisfy joint chance constraints with high probability. We present case studies that analyze the performance of the proposed approach.


Introduction
The optimization of chemical processes presents distinctive challenges to the stochastic systems community given that they suffer from three conditions: 1) there is no precise known model for most industrial scale processes (plant-model mismatch), leading to inaccurate predictions and convergence to suboptimal solutions, 2) the process is affected by disturbances (i.e. it is stochastic), and 3) state constraints must be satisfied due to operational and safety concerns, therefore constraint violation can be detrimental or even dangerous. In this work we use constrained policy search, a reinforcement learning (RL) technique, to address the above challenges.
RL is a machine learning technique that computes a policy which learns to perform a task by interacting with the (stochastic) environment. RL has been shown to be a powerful control approach, and one of the few control techniques able to handle nonlinear stochastic optimal control problems [1,2]. RL in the approximate dynamic programming (ADP) sense has been studied for chemical process control in [3]. A model-based strategy and a model-free strategy for control of nonlinear processes were proposed in [4]. ADP strategies were used to address fed-batch reactor optimization, in [5] mixed-integer decision problems were addressed with applications to scheduling. In [6] RL was combined with distributed optimization techniques, to solve an input-constrained optimal control, among other works (e.g. [7], [8]). All these approaches rely on action-value methods, which approximate the solution of the HamiltonJacobiBellman equation (HJBE), and have been shown to be reliable for some problem instances. RL also shares similar features with multiparametric model predictive control (or explicit MPC) [9,10], where an explicit representation of the controller can be found, such that it satisfies the optimality conditions. However, explicit MPC usually requires some approximations or assumptions on the models.
Policy gradient RL methods [11] have been proposed to directly optimize the control policy. Unlike action-value RL methods where convergence to local optima is not guaranteed, policy gradient methods can guarantee convergence to local optimality even in high dimensional continuous state and action spaces. However, the inclusion of constraints in policy gradient methods is not straightforward. Existing methods cannot guarantee strict feasibility of the policies even when initialized with feasible initial policies [12]. The main approaches to incorporate constraints make use of trust-region, fixed penalties [13,14], and cross entropy [12]. Furthermore, if online optimization is to be avoided, addressing the constraints by the use of penalties is a natural choice. Various approaches have been proposed in this direction, however, current approaches easily lose optimality or feasibility [13] and guarantee feasibility only in expectation. Following this thread of thought, a Lyapunov-based approach is implemented in [15], where a Lyapunov function is constructed and the unconstrained policy is projected to a safety layer allowing the satisfaction of constraints in expectation. In [16] an upper bound on the expected constraint violation is provided for the projected-based policy optimization, where the initial unconstrained policy is projected back to the constraint set. In [17] an interior-point inspired method that is widely used in control [18] is proposed, where constraints are incorporated into the reward using a logarithm barrier function, allowing the satisfaction of the constraints in expectation.
The above methods all guarantee constraint satisfaction in expectation, which is inadequate for safety critical engineering applications (very loosely speaking, this means violating 50% of the time). Furthermore, although penalties are a natural way to avoid an online optimization loop (one of the main advantages of policy gradients), tuning these penalties is not always straightforward, and usually rely on heuristics which significantly affect the performance of policy gradient algorithms [19]. This suggests constrained policy gradient methods would benefit from self-tuning parameters. As mentioned earlier, it is well known that policy gradient methods present many advantages, however, their application domain will remain limited until they can handle constraint satisfaction reliably. This is the main challenge addressed in this work. Our proposed method -chance constrained policy optimization (CCPO) -guarantees the satisfaction of joint chance constraints for the optimal policy. This allows for the satisfaction of constraints with a high probability, rather than only in expectation. To achieve the satisfaction of joint chance constraints without the need of an online optimization, we use tightening of constraints, which are appended to the objective. We treat the problem of finding the least tightening that satisfies the probabilistic constraints as a black-box optimization problem, and can therefore solve it efficiently by existing methods. Upon convergence the proposed algorithm finds an optimal policy which guarantees the constraint satisfaction to the desired tolerance.
The structure of this paper is as follows. The problem statement is outlined in section 2, the details of the proposed method for probabilistic satisfaction in RL is presented in section 3. A case study is presented in section 4, where the framework is applied to a dynamic bioreactor system, and in the last section, conclusions are outlined.

Problem Statement
In this work, the dynamic system is assumed to be given by a probability distribution, following a Markov process, where p(x i ) is the probability density function of x i , with x ∈ R nx representing the states, u ∈ R nu the control inputs and t discrete time. This behaviour is observed in systems when stochastic disturbances are present and/or other uncertainties affect the physical system, like parametric uncertainties. A discrete-time system with disturbances and parametric uncertainties can be written as: where w ∈ R nw is a vector of disturbances and p ∈ R np are uncertain parameters. This can be represented by (1). In this work we seek to maximize an objective function in expectation by using an optimal stochastic policy subject to probabilistic constraints despite the uncertainty of the system. This problem can be written as a stochastic optimal control problem (SOCP): where J is the objective function, U denotes the set of hard constraints for the control inputs, X i represents constraints for states that must be satisfied with a probability 1 − α. Specifically, with g j,t being the j th constraint to be satisfied at time instant t and the joint chance (also called probabilistic) constraints (P( T t=0 {x t ∈ X i }) ≥ 1 − α) are satisfied for the full trajectory over all t ∈ {0, . . . , T }. The probability density function of u t given state x t is p(u t |x t ), and π(x t ) is the state feedback policy. Additionally, π(·) is the stochastic policy. Unfortunately, this SOCP is generally intractable and approximations must be sought, hence the use of RL [1,20,21].
In RL a policy π θ (·) parametrized by the parameters θ is constructed. This policy maximizes the expectation of the objective function J(·). In the finite horizon discrete-time case, this objective function can be defined as where G ∈ [0, 1] is the discount factor and R t a given reward at the time instance t for the values of u t , x t .
Notice that we seek for the policy to satisfy joint chance constraints with some high probability. Previous approaches have address the satisfaction of constraints in expectation; this means that constraints are violated (roughly speaking) half of the time.

Constrained Policy Optimization for Chance Constraints
In RL agents take actions to maximize some expected reward given a performance metric. In process control, these agents become the controller, which use a feedback policy π(·), to optimize the expected value of an economic metric of the process (J). The physical system (or the model) at each sampling time produces a value for the reward R which reflects the performance of the policy. The RL algorithm determines the feedback policy that produces the greatest reward in expectation, which is referred to as policy optimization. However, there is no natural way to handle constraints, and a constraint satisfied in expectation may still have a very high probability of not being satisfied. To satisfy the constraints with some high probability and not only in expectation, we tighten the constraint with backoffs [22,23] b j,t as: where variables b j,t ≥ 0 represent the backoffs which tighten the original feasible set X t defined in (4). Backoffs restrict the perceived feasible space by the controller, and allow guarantees on the satisfaction of chance constraints. We denote τ τ τ , as the joint random variable of states, controls and rewards for a trajectory with a time horizon T : We also assume that the policy can be parametrized by a finite number of parameters θ, and denote this parametrized policy as: where u t ∈ U is inherently satisfied by the construction of the policy, e.g. the policy passes through a bounded and differential squashing function [24], and D t is a window of past inputs and states that are used by the policy, for example, if the parametrized policy is a recurrent neural network, then D t corresponds to a number of past states and controls. Recently a methodology was proposed that satisfies the expected value of the constraints [14,12], however this is not adequate for safety critical constraints in chemical processes. Instead, to account for constraint violations, we incorporate probabilistic constraints. The problem can then be reformulated as: where p(τ τ τ |θ) represents the probability of the trajectory τ τ τ given the parametrization θ of the feedback policy, and π θ * is the optimal policy.
In order to solve (9) we propose to: 1. Parameterize the feedback policy by a multilayer neural network that computes the mean and variance of the control actions, which are consequently sampled as a normal distribution resulting in a stochastic policy.
2. The probabilistic constraint in (9) is substituted by a tightened constraint set X t to guarantee closed-loop probabilistic constraint satisfaction.
3. The tightened constraints g j,t (x t ) + b j,t ≤ 0, j = 1, . . . , n g are incorporated into the objective function and can be incorporated in previous approached [13,25]. This avoids the need to explicitly solve an optimization problem at run-time. 4. The policy optimization is performed by a policy gradient framework [11].
5. The value of the backoffs should be the smallest value that guarantees constraint satisfaction with a probability of at least 1 − α. This is formulated as a black-box optimization problem, which can be solved efficiently by existing methods [26,27,28].
Notice that the value of the backoffs imply a trade-off: large values guarantee constraint satisfaction, but they make the problem over-conservative and mitigate performance, while smaller values produce solutions with high rewards, but may not guarantee the constraint satisfaction to the desired accuracy (1 − α). Also notice that if backoffs are too large, the problem might become infeasible.
In the next subsections we introduce the components that were outlined above.

Policy parametrization
For the policy parametrization we use recurrent neural networks, (RNNs) [29,30]. Let N be the length of the window of past states and controls to be used by the parametrized feedback policy, i.e.
Then the stochastic policy can be defined as: where s θ is the multilayer RNN parametrized by θ, and µ µ µ u t and Σ Σ Σ u t are the mean and covariance of the normal distribution from which the control u t is drawn. Deep structures are employed to enhance the performance of the learning process [31,32].

Probabilistic constraints
In this section, we present a strategy to handle (joint) chance constraints by policy gradient methods. There is no closed form solution for general probabilistic constraints, and therefore we compute their empirical cumulative distribution function (ECDF) instead. Furthermore, because ECDFs are themselves random variables, we must guarantee their satisfaction with some confidence, in this case we wish to define a high confidence 1 − . To guarantee with high probability (and not only in expectation) the satisfaction of constraints, we introduce constraint tightening by using backoffs b j,t (see Eq. (6)). Therefore, we conduct closed-loop Monte Carlo (MC) simulations to compute the backoffs. We pose the problem of attaining the smallest backoffs (least restrictive) that still satisfy our constraints with a probability of at least 1 − α as an expensive black-box optimization problem, and is solved by Bayesian optimization. In this way we compute a policy that satisfies our constraints with a probability of at least 1 − α, and we can be certain of this with a confidence of at least 1 − . The rest of this subsection details the aforementioned methodology.
For convenience we define a single-variate random variable C(·) representing the satisfaction of the joint chance constraint: where X = [x 1 , . . . , x T ] T , and F is the cumulative distribution function (CDF). There is no analytical expression for general probabilistic constraints, and we therefore approximate them by a non-parametric approximation, i.e. their empirical cumulative distribution function (ECDF). We therefore compute the sample approximation of the ECDF by S Monte Carlo (MC) simulations of X: where X s = [x s 1 , . . . , x s T ] T is the s th MC sample of the state trajectory,F S is the approximate probability for a trajectory to satisfy all constraints, and 1{C(X) ≤ 0} is the indicator function for a single trajectory to satisfy all constraints defined as follows, The indicator function is a Bernoulli random variable, which means thatF S follows a Binomial distribution, withF S ∼ 1 S Bin(S, F ), F being the cumulative distribution function (CDF). The confidence bound for the ECDF can then be computed from the Binomial CDF [33].
Notice that the quality of the approximation in (13) strongly depends on the number of samples S used and it is therefore desirable to quantify the uncertainty of the sample approximation itself. This problem has been studied to a great extent in the statistics literature [33], leading to the following Lemma 1.

Lemma 1 ([33]
). Assume we are given a value of the ECDFF S (see (13)) based on S i.i.d. samples, then the true value of the CDF, F , has the following lower boundF lb with a confidence level of 1 − : with betainv(·, ·, ·) being the inverse of the beta CDF with parameters {S + 1 − SF S } and {SF S }.
Lemma 1 states that it is possible to compute a probabilistic lower bound, with an arbitrarily confidence of at least 1 − , for the ECDF of the chance constraints. Using the above result, we can state the following theorem.
Theorem 2. Given a value of the ECDFF S based on S independent samples and a corresponding lower boundF lb ≥ 1 − α. Then, with a confidence level of 1 − , our original chance constraint optimization problem in Eq. (9) holds.
Proof. From S sampled trajectories of X a valid lower boundF lb to the true CDF value F can be determined from Theorem 1 with a confidence level of 1 − . IfF lb is greater than or equal to 1 − α, then the following probabilistic bound holds on the true CDF value F according to Theorem The implication of Theorem 2 is that we can guarantee the satisfaction of joint chance constraints with a user-defined probability of at least 1 − α and a user-defined confidence of at least 1 − . Based on Theorem 2 we now wish to obtain a value forF lb so that the original constraints (P( are satisfied with a confidence of 1 − . LargerF lb lead to more conservative solutions, and consequently worse values for the objective function. Therefore, the least conservative solution leading to the same guarantees is realized ifF lb = 1 − α given thatF lb is the tightest lower bound for F . Hence, we aim to compute tightened constraints employing backoffs b j,t such thatF lb = 1 − α is satisfied when the policy optimization has terminated. Theorem 2 allows us to find a solution to the original problem having satisfied all joint chance constraints with a confidence at least 1 − . To compute the tightened constraint set as denoted in Eq. (6), we first compute an initial set of backoffs (b 0 j,t ), as it has been proposed in [34], where with δ being a tuning parameter for the initial backoffs b 0 j,t . Additionally, the expected value of g j,t is approximated with a sample average approximation (SAA): The initial backoff values are computed to probabilistically satisfy each constraint: with g j,t = [g j,t (x 1 t ), . . . , g j,t (x S t )] T and F −1 gj,t (1 − δ) is the 1 − δ quantile of g j,t (17a). We wish the backoffs b (with elements b j,t ∀ j, t ∈ {1, . . . , n g } × {1, . . . , T }) to be large enough to ensure the probabilistic satisfaction of constraints. However, if the backoffs are too large it may result in a conservative solution, therefore a worse performance, or even infeasibility of the problem. To obtain the least conservative solution that still guarantees the probabilistic constraint satisfaction specified by 1 − α we solve a root-finding problem using Theorem 1 to find a vector γ γ γ that parametrizes b j,t : Notice, that nowF lb is a function of γ γ γ = [γ 1 , ..., γ ng ]. With the above procedure we compute a deterministic surrogate for the constraints that allows us to satisfy the joint chance constraints in the original optimization problem (9). In the subsequent section we explain how this constraint surrogate is incorporated into the reinforcement learning framework.

Policy gradient for fixed backoffs
In this work we reformulate chance constraints such that their satisfaction is guaranteed with high probability, and they can be incorporated into the objective of the policy gradient method [11]. We therefore avoid the need for a numerical optimization every time the agent/controller outputs an action/control input.
Loosely speaking policy gradient methods aim to update the parametrized policy using the gradient of the reward. Without loss of generality, we outline the implementation for the REINFORCE [35] algorithm for ease of presentation, but this can be generalized to any policy gradient or actor-critic method [20,36].
Previous works have embedded Lagrangian methods as adaptive penalty coefficients to enforce satisfaction of constraints [14], however, an adaptive scheme such as gradient ascent or its variants on inequality multipliers are difficult to justify in theory (see for example [37] chapter 12 or [38] chapter 5), and in practice tend to have numerical issue when an arbitrary number of constraints are enforced and different constraints are active in different instances [37].
In this work, we instead propose a p-norm : where [g j,t (x t ) + b j,t ] − = max{g j,t , (x t ) + b j,t , 0}, || · || p is a p-norm of the vector g t = g 1,t , . . . , g ng,t . We advocate for the use of l 1 (p = 1) and l 2 (p = 2) norms due to the arbitrary number of constraints that may be considered and their numerical stability [37].
Theorem 3. (Probabilistic satisfaction of constraints) Consider the chance constrained stochastic optimal problem (9) and let π θ * be the trained policy that maximizes the expected reward (see [13,14] or (19)) and satisfies (18) using backoffs b. Then the joint chance constraints (11) will be satisfied with probability 1 − α and a confidence level of 1 − .
Remark 2. The parameter κ can be updated using standard techniques from constrained optimization [37].
We now use the policy gradient theorem [11] to obtain an explicit gradient estimate of the reward with respect to the parameters of our policy: where denotes the realization of the s th trajectory with R s being the reward for sample s and denotes the past states and controls used by the policy for sample s. The variance of the gradient estimate can be reduced with the aid of an action-independent baselineβ S , which does not introduce a bias [39]. A simple but effective baseline is the expectation of reward under the current policy, approximated by the mean of the sampled paths: which leads to: Using the above gradient of the expected reward with respect to the parameters the policy can now be iterative adjusted, in a steepest ascent framework or any of its variants (e.g. Adam): where k is the adaptive learning rate. The algorithm that trains the policy network for a fixed backoff value b is outlined in Algorithm 1.

Backoff iterations
Solving Eq. (18) is complicated because there is no closed form solution, or even an explicit expression. To solve the root-finding problem in Eq. (18) efficiently, we formulate it as a least-squares expensive black box optimization problem [27] and solve it via Bayesian optimization with the objective function Algorithm 1 yields an optimal stochastic policy for fixed values of the backoffs. We propose Algorithm 2 to iteratively adjust the backoffs to guarantee probabilistic constraint satisfaction. A description of the steps conducted in Algorithm 2 is presented here: Step (1): The policy is trained with b = 0 using Algorithm 1.
Step (2): The initial estimates of the backoff are computed by Monte-Carlo closed-loop simulations.
Step (4): For each m th iteration: (4i) A Gaussian process which maps the backoff values to the squared residuals in Eq. (25) is constructed. We compute the objective function value as in Eq. (25) We conduct a Bayesian optimization step, via lower confidence bound minimization [40] of the GP. This allows us to obtain the next value for γ γ γ NΓ+m .
(4iv) We compute the new backoff b with the new value of the backoff parameter γ γ γ NΓ+m . We also compute the new value for the residuals (objective) F. Data matrices Γ Γ Γ and F are updated. The algorithm goes back to (i) where a new GP is constructed incorporating the new values of γ γ γ NΓ+m and F and the algorithm proceeds until some tolerance is achieved (e.g. F(γ γ γ NΓ+m ) ≤ tol).
It should be noted that every time a new backoff is computed the policy is re-optimized. This may look inefficient at first glance, however the convergence is achieved fast, as at every iteration, the previous policy is used as an initial guess for the next iteration, making the first iteration the most expensive one. Additionally, because the problem is treated as an expensive black-box optimization problem the number of iterations from this outer loop is in general quite small (in our examples no more than 13 iterations). (see section 4.3). By the end of Algorithm 2, a probabilistically constrained policy will have been constructed.

Policy initialization
Reinforcement learning methods (particularly policy gradient) are computationally expensive; mainly because initially the agent (or controller in our case) explores the control action space randomly. In the case of process optimization and control, it is possible to use a preliminary controller, along with supervised learning to hot-start the policy, and significantly speed-up convergence. The initial parameterization for the policy (before Step (1)) is trained in a supervised learning fashion where the states are the inputs and the control actions are the outputs. This has been further discussed in [41,42,43]

Case Studies
The case studies in this paper focuses on the photo-production of phycocyanin synthesized by cyanobacterium Arthrospira platensis. The two case studies are separated regarding the type of uncertainty: The case study 1 considers parametric uncertainty and model is considered to be known. The second case study considers that only data is available, where additive disturbance and measurement noise are present (no knowledge of the system's equations). Additionally, in the two case studies different penalizations of the constraints are implemented: Case study 1 uses Eq. (19) with arbitrary κ = 0.1 and p = 1, and Case study 2 uses Eq. (19) with arbitrary κ = 1. and p = 2.
Phycocyanin is a high-value bioproduct and its biological function is to enhance the photosynthetic efficiency of cyanobacteria and red algae. It has applications as a natural colorant to replace other toxic synthetic pigments in both food and cosmetic production. Additionally, the pharmaceutical industry considers it beneficial because of its unique antioxidant, neuroprotective, and anti-inflammatory properties.
The dynamic system consists of the following system of ODEs describing the evolution of the concentration (c) of biomass (x), nitrate (N ), and product (q). The dynamic model is based on Monod kinetics, which describes microorganism growth in nutrient sufficient cultures, where intracellular nutrient concentration is kept constant because of the rapid replenishment. We assume a fixed volume fed-batch. The manipulated variables as in the previous examples are the light intensity (u 1 = I) and inflow rate (u 2 = F N ). The mass balance equations are The parameter values are adopted from [22].

Case Study 1
Uncertainty is assumed for the initial concentration, where c x (0) c N (0) ∼ N ( 1. 150. , diag(1× 10 −3 , 22.5)) and c q (0) = 0. Additionally, 10% of parametric uncertainty for the system is assumed: type of uncertainty is common in engineering settings, as the parameters are obtained using experimental data, and they are subject to their respective confidence regions after they are estimated using regression techniques. The objective function (reward) in this work is to maximize the product's concentration (c q ) at the end of the batch. The objective is additionally penalized by the change of the control actions As a result the reward is: where t ∈ {0, T − 1} and ∆u t = u t − u t−1 . The constraints in this work for each time step are c N ≤ 800 and c q ≤ 0.011c X . This constraints have been normalized as: and the joint chance constraint is meant to be satisfied with probability 99% (α = 0.01) and confidence level is 99% ( = 0.01). The constraints are added as a penalty with κ = 1 using (19) and Additionally, the Gaussian process has zero mean as prior, squared-exponential (SE) kernel as the covariance function, and the inputs-outputs are normalized using its mean and variance respectively. After the completion of the training the backoffs have been computed to satisfy (18). In a rather small number of iterations the backoff values managed to force the F lb to 0.99. Now, the actual closed-loop constraint satisfaction can be depicted in Fig. 1(a) and Fig. 2(a), where the shaded areas are the 98% and 2% percentiles. Notice, that even though the figures for both methods look similar the constraint satisfaction is significantly different, this can be observed in Fig. 1(b) and Fig. 2(b), where the region that the violation of constraints occurs has been zoomed in. This result is also quantitatively depicted in Table 1. The column labelled 'actual' is the probability of constraint violation that corresponds to the fraction of the 1000 Monte-Carlo trajectories that satisfied both of the constraints (see Eq. (13)). The column labelled 'desired' is the goal for the probability for constraint satisfaction. It is clear that when backoffs are not applied, almost 50% of the constraints are violated. On the other hand when backoffs are applied then all the constraints are satisfied, which is an expected result as the goal was the lower bound of ECDF (F lb ) to be 0.99, which means that the actual probability is equal or higher.  The backoff values for each update are shown in Fig. 3 (a, b), where the red-dashed represents the converged final value. It should be noted that the final value for the product c q (29) is 0.163 and 0.167 when backoffs are applied and when they are not. This difference is to be expected given that in the nominal case (b = 0), the policy allows the violation of constraints which lead to an increase in the concentration of the product.

Case Study 2
The second case study considers the same system but this time no equation is assumed to be available. In this case, a Gaussian process [44] is used to model the available data. Different kernels can have different effects in the performance, in this case the Matérn32 [44] kernel is employed. The data was generated using the system in Eqs. (26)(27)(28), and an additive normally distributed disturbance (w) and measurement noise (v) with zero mean and Σ Σ Σ w = diag(4 × 10 −4 , 0.1, 1 × 10 −8 ), Σ Σ Σ v = diag(4 × 10 −5 , 0.01, 1 × 10 −9 ). The same uncertainty for the initial conditions is assumed.
Bayesian frameworks have been used in reinforcement learning [45,46] as they can model both the epistemic and aleatoric uncertainty (in Gaussian processes, usually homoscedastic aleatoric uncertainty is considered). One of the fundamental steps here is the propagation of the uncertainty for each step. The nature of Gaussian processes has been exploited in [47,48,22], where the generated samples from the GP are conditioned in the posterior (without noise). Then, the GP predicts a mean µ µ µ x and variance Σ Σ Σ x and the state x is sampled from a normal distribution: x ∼ N (µ µ µ x , Σ Σ Σ x ). In this setting 8 episodes were randomly generated using a Sobol sequences [49] for the control variables of all 12 time intervals and the initial conditions to train the Gaussian process.
After the training the actual closed-loop constraint satisfaction can be validated. Fig. 4(a) and Fig. 5(a) illustrate the path constraints for each time interval, where the shaded areas are the 98% and 2% percentiles. The results are compared with the case in absence of backoffs (b = 0). It should be noticed that in the absence of backoffs the mean of the constraints' samples barely satisfies the bounds. The difference is clear and more apparent in Fig. 5(a). This can further be noticed in Fig. 4(b) and Fig. 5(b), where focus has been put on the region of constraint violation. Table 1 shows the probability of constraint violation that corresponds to the fraction of the 1000 Monte-Carlo trajectories that satisfied the constraints (see Eq. (13)) compared to the 'desired' probability for constraint satisfaction. The absence of constraint tightening results in only 24% of constraint satisfaction compare to our proposed method where 97% are satisfied. Notice that the desired probability designed in terms of the lower bound of the ECDF (F lb ) is set to be 0.95.  The backoff values for each update are shown in Fig. 6 (a, b), where the red-dashed represents the converged final value. It should be noted that the final value for the product c q (29) is 0.153 and 0.171 when backoffs are applied and when they are not.  The algorithm is implemented in Pytorch [50] version 0.4.1. Adam [51] is employed to compute the network's parameter values using a step size of 10 −2 with the rest of hyperparameters at their default values.

Conclusions
In this paper we address the problem of finding a policy that can satisfy constraints with high probability. The proposed algorithm -chance constrained policy optimization (CCPO) -uses constraint tightening by applying backoffs to the original feasible set. Backoffs restrict the perceived feasible space by the controller, and allow guarantees on the satisfaction of chance constraints. We find the smallest backoffs (least conservative) that still guarantee the desired probability of satisfaction by stating the root-finding problem as a black-box optimization problem. This allows the algorithm to construct a policy that can guarantee the satisfaction of joint chance constraints with a user-defined probability of at least 1 − α and a user-defined confidence of at least 1 − . Furthermore, the proposed methodology can be combined with other penalty or Lagrangian approaches for constrained policy search. Being able to solve constraint policy optimization problems with high probability satisfaction has been one of the main bottlenecks for the wider use of reinforcement learning in engineering applications. This work aims to take a step towards applying RL to the real world, where constraints on policies are necessary for the sake of safety and product quality.

Appendix A. Bayesian optimization
Bayesian optimization is a popular approach for black-box optimization [27]. A review of Bayesian optimization can be found in [52].
In this paper Bayesian optimization is exploited to obtain the minimizer γ γ γ * of F(γ γ γ) = F lb − (1 − α) 2 : γ γ γ * ∈ arg min γ γ γ F(γ γ γ) (A.1) The function F(γ γ γ) to be minimized can only be observed through an unbiased noisy observation based on MC estimates ofF lb : where ω is assumed to be Gaussian white noise. The noise ω is assumed to be unknown and estimated within the GP framework. Commonly Gaussian processes (GP) are employed as nonparametric models.
To start the algorithm we first require some data to build the initial GP, i.e. the function F(·) is queried at N Γ points. In general the input data-points Γ Γ Γ = [γ γ γ 1 , . . . , γ γ γ NΓ ] are selected based on a space-filling design, such as a Latin hybercube design [53]. From this we obtain the corresponding responsesF lb = [F(γ γ γ 1 ), . . . , F(γ γ γ NΓ )]. A GP model can then be trained from the input-output data. In particular the hyperparameters of the GP were trained using maximum likelihood estimation. The GP model can then be utilized to obtain the Gaussian distribution of F(γ γ γ) at an arbitrary query point γ γ γ [44]: where µ GP (γ γ γ) and σ 2 GP (γ γ γ) are the mean and variance prediction of the GP respectively. In this setting the approach sequentially selects a location γ γ γ at which to query F(·) and observe y. After a selected number of iterations M the algorithm returns a best-estimate of γ γ γ * . To accomplish this the GP of F(·) is iteratively updated from the available data of F(·). One could simply sample at the minimum of the mean function µ GP (γ γ γ). However, sampling at a point with higher uncertainty could yield a lower minimum, i.e. there is a trade-off between sampling at points with low values of the mean function µ GP (γ γ γ) and high values of the variance function σ 2 GP (γ γ γ). The selection of the query points is given by so-called acquisition function, for which we used the lower confidence bound: γ γ γ m = arg min γ γ γ µ GP (γ γ γ) − 3σ GP (γ γ γ) (A. 4) where γ γ γ m denotes the query point at iteration m. Note that the query point at each iteration are chosen at points that are predicted to be low by the mean function, but could also potentially yield lower values according to the variance function.