Statistical Learning for Probability-Constrained Stochastic Optimal Control

TL;DR

This paper develops Monte Carlo algorithms using statistical learning methods like logistic and Gaussian process regression to solve stochastic control problems with probabilistic constraints, demonstrated on microgrid management.

Contribution

It introduces a Regression Monte Carlo framework incorporating statistical tools for learning admissible control sets under probabilistic constraints, improving efficiency and reliability.

Findings

Logistic and Gaussian process regression outperform other methods in estimating admissibility probabilities.

The proposed algorithms extend RMC to handle probability constraints effectively.

Case studies validate the approach in microgrid control scenarios.

Abstract

We investigate Monte Carlo based algorithms for solving stochastic control problems with probabilistic constraints. Our motivation comes from microgrid management, where the controller tries to optimally dispatch a diesel generator while maintaining low probability of blackouts. The key question we investigate are empirical simulation procedures for learning the admissible control set that is specified implicitly through a probability constraint on the system state. We propose a variety of relevant statistical tools including logistic regression, Gaussian process regression, quantile regression and support vector machines, which we then incorporate into an overall Regression Monte Carlo (RMC) framework for approximate dynamic programming. Our results indicate that using logistic or Gaussian process regression to estimate the admissibility probability outperforms the other options. Our algorithms offer an efficient and reliable extension of RMC to probability-constrained control. We illustrate our findings with two case studies for the microgrid problem.

Figures

Tables

Table 1. Table 1 : Parameters for the Microgrid example.

$I_{\max }=10$ (kWh), $B_{\min }=-6,B_{\max }=6$ (kW), $𝒦=5$

$T=48$ (hours), $\Delta t=0.25$ (hours)

Table 2. Table 2 : Parameters for the estimation of the admissible sets for each method. We use total simulation budget of 10 5 superscript 10 5 10^{5} for all models except the Gold Standard.

Method	Budget ($M_a\times M_b$)	Further parameters
Gaussian Process (GPR)	$2000\times 50$	Matern-3/2 kernel
Logistic Regression (LR)	${10}^5\times 1$	Degree-2 polynomials
Parametric Density Fitting (PF)	$2000\times 50$	Truncated Gaussian, Matern-3/2 kernel
Empirical Percentile (EP)	$1000\times 100$	Squared exponential kernel
Conditional Tail Expectation (CTE)	$1000\times 100$	Squared exponential kernel
Quantile Regression (QR)	${10}^5\times 1$	Degree-4 polynomials
Support Vector Machine (SVM)	$2000\times 50$	C =1, RBF kernel
Gold Standard (GS)	${10}^6\times {10}^4$	budget = ${10}^{10}$

Table 3. Table 3 : Cost of running the microgrid V ^ 0 ​ ( 0 , 5 ) subscript ^ 𝑉 0 0 5 \hat{V}_{0}(0,5) , frequency of inadmissible decisions w f ​ r ​ e ​ q subscript 𝑤 𝑓 𝑟 𝑒 𝑞 w_{freq} , average inadmissibility margin w a ​ v ​ m subscript 𝑤 𝑎 𝑣 𝑚 w_{avm} , realized frequency of violations (i.e. blackouts) w r ​ l ​ z ​ d subscript 𝑤 𝑟 𝑙 𝑧 𝑑 w_{rlzd} , test statistic 𝒯 ~ ~ 𝒯 \tilde{\mathcal{T}} and frequency of binding constraint w b ​ i ​ n ​ d subscript 𝑤 𝑏 𝑖 𝑛 𝑑 w_{bind} for the example in Section  5.2 .

Method	${\hat{V}}_0(0,5)$ ($)	$w_{freq}$ (%)	$w_{avm}$ (kW)	$w_{rlzd}$ (%)	$\widetilde{𝒯}$	$w_{bind}$ (%)
GS	26.79	0.00	0.00	0.37	-	-
LR	26.83	0.09	0.82	0.03	-125	8.69
GPR	26.89	0.53	0.16	0.11	-98	8.10
PF	26.79	1.36	0.27	0.21	-69	8.51
SVM	26.68	5.26	0.55	1.83	388	9.67
QR	27.04	5.95	0.33	0.98	145	9.49
CTE	26.99	7.79	0.43	1.63	320	9.93
EP	26.36	8.39	0.49	1.98	403	10.45

Table 4. Table 4 : Impact of margin of error ξ 𝜉 \xi on the estimated cost of running the microgrid V ^ 0 ​ ( 0 , 5 ) subscript ^ 𝑉 0 0 5 \hat{V}_{0}(0,5) , frequency of inadmissible decisions w f ​ r ​ e ​ q subscript 𝑤 𝑓 𝑟 𝑒 𝑞 w_{freq} , and test statistic 𝒯 ~ ~ 𝒯 \tilde{\mathcal{T}} from ( 31 ). The probabilistic constraint is p = 5 % 𝑝 percent 5 p=5\% .

	$\xi =0\%$			${\xi }^{\left(0.95\right)}(𝐱,u)$			$\xi =4\%$
Method	${\hat{V}}_0(0,5)$	$w_{freq}$	$\widetilde{𝒯}$	${\hat{V}}_0(0,5)$	$w_{freq}$	$\widetilde{𝒯}$	${\hat{V}}_0(0,5)$	$w_{freq}$	$\widetilde{𝒯}$
GS	26.79	0.00	-	-	-	-	-	-	-
LR	26.83	0.09	-125	26.95	0.08	-124	27.86	0.04	-112
GPR	26.89	0.53	-98	28.00	0.01	-110	28.12	0.00	-107
PF	26.79	1.36	-69	-	-	-	27.91	0.44	-96
SVM	26.68	5.26	388	29.65	3.41	225	29.60	3.41	225
QR	27.04	5.95	145	26.89	5.17	72	28.61	0.00	-117
CTE	26.99	7.79	320	27.36	7.52	274	28.44	6.83	248
ER	26.36	8.39	403	26.97	7.78	225	28.13	7.08	283

Table 5. Table 5 : Cost of running the microgrid V ^ 0 ​ ( μ ​ ( 0 ) , 5 ) subscript ^ 𝑉 0 𝜇 0 5 \hat{V}_{0}(\mu(0),5) , frequency of inadmissible decisions w f ​ r ​ e ​ q subscript 𝑤 𝑓 𝑟 𝑒 𝑞 w_{freq} , average inadmissibility margin w a ​ v ​ m subscript 𝑤 𝑎 𝑣 𝑚 w_{avm} , realized frequency of violations w r ​ l ​ z ​ d subscript 𝑤 𝑟 𝑙 𝑧 𝑑 w_{rlzd} and frequency of the constraint being binding w b ​ i ​ n ​ d subscript 𝑤 𝑏 𝑖 𝑛 𝑑 w_{bind} for the case study in Section  5.3 .

Method	${\hat{V}}_0(\mu \left(0\right),5)$	$w_{freq}$ (%)	$w_{avm}$ (kW)	$w_{rlzd}$ (%)	$\widetilde{𝒯}$	$w_{bind}$ (%)
GS	53.38	0	0	0.30	-	-
LR	53.78	0.03	0.79	0.01	-301	45.2
GPR	54.04	1.17	0.14	0.19	-220	31.0
PF	54.55	0.02	0.26	0.01	-226	25.7
SVM	40.52	43.37	0.91	43.37	5,306	46.4
QR	52.56	42.87	0.28	38.41	4,772	46.3
CTE	53.02	21.62	0.21	10.43	1,079	46.0
EP	52.82	21.91	0.23	11.57	1,227	46.1

Table 6. Table 6 : Impact of simulation budget on performance of SVM for the case study in Section  5.2 and probability thresholds p = 5 % 𝑝 percent 5 p=5\% and p = 1 % 𝑝 percent 1 p=1\% . The reported values are averages over 10 runs of each scheme. The total simulation budget is divided into batch size M b subscript 𝑀 𝑏 M_{b} and number of design sites M a subscript 𝑀 𝑎 M_{a} . For total budget 10 5 superscript 10 5 10^{5} : ( M b , M a ) = ( 100 , 1000 ) subscript 𝑀 𝑏 subscript 𝑀 𝑎 100 1000 (M_{b},M_{a})=(100,1000) ; for 10 6 superscript 10 6 10^{6} : ( M b , M a ) = ( 500 , 2000 ) subscript 𝑀 𝑏 subscript 𝑀 𝑎 500 2000 (M_{b},M_{a})=(500,2000) ; for 10 7 superscript 10 7 10^{7} : ( M b , M a ) = ( 2000 , 5000 ) subscript 𝑀 𝑏 subscript 𝑀 𝑎 2000 5000 (M_{b},M_{a})=(2000,5000) ; for 10 8 superscript 10 8 10^{8} : ( M b , M a ) = ( 10000 , 10000 ) subscript 𝑀 𝑏 subscript 𝑀 𝑎 10000 10000 (M_{b},M_{a})=(10000,10000) .

$p$	Budget	${\hat{V}}_0(0,5)$ ($)	$w_{freq}$ (%)	$w_{avm}$ (kW)	$w_{rlzd}$ (%)	$\widetilde{𝒯}$	$w_{bind}$ (%)
5%	${10}^5$	26.38	5.93	0.78	2.80	665	9.73
	${10}^6$	26.55	5.28	0.55	1.84	386	9.77
	${10}^7$	26.68	4.96	0.53	1.64	330	9.75
	${10}^8$	26.79	1.50	0.27	0.30	-51	9.22
1%	${10}^5$	28.32	6.63	0.93	2.43	1,460	9.87
	${10}^6$	28.26	5.17	0.66	1.09	631	9.56
	${10}^7$	28.52	0.55	0.24	0.03	-39	8.78
	${10}^8$	28.41	0.15	0.22	0.01	-51	8.82

Table 7. Table 7 : Impact of conservative 𝒰 ( ρ ) superscript 𝒰 𝜌 \mathcal{U}^{(\rho)} estimators for the case study in Section  5.2 . The probabilistic constraint is set at p = 5 % 𝑝 percent 5 p=5\% .

	$\rho =90\%$			$\rho =99\%$			$\rho =99.95\%$
Method	${\hat{V}}_0(0,5)$	$w_{freq}$	$w_{rlzd}$	${\hat{V}}_0(0,5)$	$w_{freq}$	$w_{rlzd}$	${\hat{V}}_0(0,5)$	$w_{freq}$	$w_{rlzd}$
LR	26.74	0.090	0.034	26.87	0.085	0.032	27.04	0.085	0.026
GPR	27.34	0.012	0.055	28.06	0.007	0.037	28.06	0.005	0.029
SVM	27.35	4.975	1.732	29.20	3.481	1.117	29.72	3.395	1.088
QR	27.20	5.373	0.793	27.18	4.880	0.676	27.04	4.409	0.591
CTE	27.93	7.158	1.153	28.31	6.766	0.888	28.61	6.163	0.714
EP	26.78	7.990	1.497	27.17	7.629	1.183	27.96	7.102	0.956

Equations

$$\mathrm{d}\mathbf{X}(t)=b(t,\mathbf{X}(t),u(t))\mathrm{d}t+\sigma (t,\mathbf{X}(t),u(t))\mathrm{d}B(t),$$

$$\mathcal{U}_{n:N}(\mathbf{X}_{n})=\Big{\{}(u_{k})_{k=n}^{N}:P_{k}(\mathbf{X}_{k},u_{k})\in\mathcal{A}_{k}\ \forall k\in\{n,\,\dots,\,N-1\}\Big{\}},$$

$$P_{n}(\mathbf{X}_{n},u_{n})\equiv p_{n}(\mathbf{X}_{n},u_{n}):=\mathbb{P}\Big{(}{\color[rgb]{0,0,0}\mathcal{G}_{n}((\mathbf{X}(s))_{s\in[t_{n},t_{n+1})})}>0\big{|}\,\mathbf{X}_{n},u_{n}\Big{)}\ \text{ and }\ \mathcal{A}_{n}:=[0,p),$$

$$\mathcal{U}_{n}(\mathbf{X}_{n}):=\Big{\{}u\in\mathcal{W}:\ p_{n}(\mathbf{X}_{n},u_{n})<p\Big{\}}.$$

$$V_{n}(\mathbf{X}_{n})=\inf_{(u_{s})_{s=n}^{N}\in\mathcal{U}_{n:N}(\mathbf{X}_{n})}\Big{\{}\mathbb{E}\Big{[}\sum_{k=n}^{N-1}\int_{t_{k}}^{t_{k+1}}\pi_{s}(\mathbf{X}(s),u_{k})\mathrm{d}s+W(\mathbf{X}(t_{N}))\Big{|}\,\mathbf{X}_{n}\Big{]}\Big{\}},$$

$$\begin{split}V_{n}(\mathbf{X}_{n})=&\inf_{u\in\mathcal{U}_{n}(\mathbf{X}_{n})}\Big{\{}\mathcal{C}_{n}(\mathbf{X}_{n},u)\Big{\}},\\ \text{where }\quad\mathcal{C}_{n}(\mathbf{X}_{n},u)&=\mathbb{E}\left[\int_{t_{n}}^{t_{n+1}}\pi_{s}(\mathbf{X}(s),u)\mathrm{d}s+V_{n+1}(\mathbf{X}(t_{n+1}))\Big{|}\,\mathbf{X}_{n},u\right].\end{split}$$

$$G_{n}(\mathbf{X}_{n},u_{n}):=\mathcal{L}\Big{(}\mathcal{G}_{n}((\mathbf{X}(s))_{s\in[t_{n},t_{n+1})})\Big{|}\,\mathbf{X}_{n},u_{n}\Big{)},$$

$$q_{n}(\mathbf{X}_{n},u_{n}):\,(\mathbf{X}_{n},u_{n})\mapsto\text{arg}\inf_{z}\Big{\{}\mathbb{P}\Big{(}G_{n}(\mathbf{X}_{n},u_{n})>z\Big{)}\leq p\Big{\}}.$$

$${\mathcal{U}}_n({\mathbf{X}}_n):=\left\{u:p_n({\mathbf{X}}_n,u)<p\right\}=\left\{u:q_n({\mathbf{X}}_n,u)\le 0\right\},$$

$$u^{\min}_{n}(\mathbf{X}_{n}):=\inf_{u\in\mathcal{W}}\Big{\{}u:p_{n}(\mathbf{X}_{n},u)<p\Big{\}}.$$

$$\Big{\{}u\in\mathcal{W}:\mathbb{E}\Big{[}g\Big{(}G_{n}(\mathbf{X}_{n},u)\Big{)}\Big{]}\leq p\Big{\}},$$

$$\mathcal{X}^{a}_{n}(u):=\Big{\{}\mathbf{x}\in\mathcal{X}:p_{n}(\mathbf{x},u)<p\Big{\}}.$$

$$\begin{split}\hat{V}_{n}(\mathbf{X}_{n})=&\inf_{u_{n}\in\hat{\mathcal{U}}_{n}(\mathbf{X}_{n})}\Big{\{}\hat{\mathcal{C}}_{n}(\mathbf{X}_{n},u_{n})\Big{\}},\\ \text{where }\quad\hat{\mathcal{C}}_{n}(\mathbf{X}_{n},u_{n})&:=\hat{\mathbb{E}}\left[\int_{t_{n}}^{t_{n+1}}\pi_{s}(\mathbf{X}(s),u_{n})\mathrm{d}s+\hat{V}_{n+1}(\mathbf{X}(t_{n+1}))\Big{|}\mathbf{X}_{n},u_{n}\right].\end{split}$$

$$\psi\big{(}(\mathbf{X}(s))_{s\in[t_{n},t_{n+1}]}\big{)}=\int_{t_{n}}^{t_{n+1}}\pi_{s}(\mathbf{X}(s),u_{n})\mathrm{d}s+\hat{V}_{n+1}(\mathbf{X}(t_{n+1})).$$

$${\bar{p}}_n(\mathbf{x},u):=\sum _{b=1}^{M_b}\frac{\text{mathbbm}{1}_{g_n^b(\mathbf{x},u)>0}}{M_b}$$

$$\mathcal{U}_{n}(L_{n},I_{n},C_{n}):=\left\{u:\mathbb{P}\Big{(}\sup_{s\in[t_{n},t_{n+1})}S(s)>0\Big{|}(L_{n},I_{n},C_{n},u)\Big{)}<p\right\}.$$

$$G_n(L_n,I_n,C_n,u_n)$$

$$p_n(L_n,I_n,C_n,u_n)$$

$$\displaystyle V_{n}(L_{n},I_{n},C_{n})=\underset{\{u_{k}\}_{k=n}^{N-1}}{\min}\left\{\mathbb{E}\left[\sum_{k=n}^{N-1}\Big{(}\mathbbm{1}_{\{{C}_{k}=0,u_{k}>0\}}\mathcal{K}+\rho(u_{k})\Delta t_{k}\Big{)}+W(L_{N},I_{N},C_{N})\Big{|}(L_{n},I_{n},C_{n})\right]\right\},$$

$$\displaystyle\text{subject to }\qquad\mathbb{P}\Big{(}\underset{s\in[t_{k},t_{k+1})}{\sup}S(s)>0\Big{|}(L_{k},I_{k},C_{k},u_{k})\Big{)}<p,\qquad\,k=n,\ldots,N-1,$$

$$\text{mathbbm}{1}_{\{C_n=0,u_n>0\}}\mathcal{K}+\rho (u_n)\Delta t_n.$$

$$\begin{array}{rl}\text{Generate design}& \to \text{Generate 1-step paths \& admissibility statistic}\to \text{Estimate admissible set}\\ \text{Generate design}& \to \text{Generate 1-step paths \& pathwise profits}\to \text{Estimate continuation function}\end{array}$$

$$w_{n}^{i}:=\mathbbm{1}\Big{(}\mathcal{G}((\mathbf{x}_{n_{k}}^{i})_{k\in\{0,\ldots,K-1\}})>0\Big{)},\qquad i=1,\ldots,M_{a},$$

$${\hat{p}}_n:=\arg \underset{f_n^a\in {\mathcal{H}}_n^a}{\min }\sum _{i=1}^{M_a}{\mathcal{L}}_n^a(f_n^a,w_n^i;{\mathbf{x}}_n^i,u_n^i).$$

$${\hat{p}}_n^{(\rho )}(\mathbf{x},u)$$

$$:={\bar{p}}_n(\mathbf{x},u)+z_{\rho }\sqrt{\frac{{\bar{p}}_n(\mathbf{x},u)(1-{\bar{p}}_n(\mathbf{x},u))}{M_b}},$$

$${\hat{\mathcal{U}}}_n^{(\rho )}(\mathbf{x}):={\hat{\mathcal{U}}}_n^{{\xi }^{(\rho )}}(\mathbf{x})=\left\{u:{\hat{p}}_n(\mathbf{x},u)+{\xi }_n^{(\rho )}(\mathbf{x},u)<p\right\}.$$

$${\hat{\mathcal{U}}}_n^{\xi }(\mathbf{x})=\left\{u:{\hat{p}}_n(\mathbf{x},u)+{\xi }_n(\mathbf{x},u)<p\right\},$$

$${\hat{\mathcal{U}}}_n^{\xi =c}(\mathbf{x})=\left\{u:{\hat{p}}_n(\mathbf{x},u)+c<p\right\}.$$

$$y_n^j=\sum _{k=0}^{K-1}{\pi }_{n_k}({\mathbf{x}}_{n_k}^j,u_n^j)\Delta n_k+v_{n+1}^j,\text{ where }{v}_{n+1}^j=\underset{u\in {\hat{\mathcal{U}}}_{n+1}({\mathbf{x}}_{n_K}^j)}{\inf }{\hat{\mathcal{C}}}_{n+1}({\mathbf{x}}_{n_K}^j,u),j=1\dots M_c,$$

$${\hat{\mathcal{C}}}_n(\cdot ,\cdot ):=\arg \underset{f_n^c\in {\mathcal{H}}_n^c}{\min }\sum _{n=1}^{M_c}|f_n^c({\mathbf{x}}_n^j,u_n^j)-y_n^j{|}^2.$$

$$\hat{\bm{\alpha}}:=\arg\min_{\bm{\alpha}}\sum_{j=1}^{M_{c}}\Big{|}\sum_{k}\alpha_{k}\phi_{k}(\mathbf{x}^{j},u^{j})-y^{j}\Big{|}^{2}.$$

Full text

{AMS}

93E20, 93E35, 49L20

1 Introduction

Stochastic control with probabilistic constraints is a natural relaxation of deterministic restrictions which tend to generate high costs forcing the avoidance of extreme events no matter their likelihood of occurrence. In contrast, with probabilistic constraints, constraint violation is tolerated up to a certain level offering a better trade-off between admissibility and cost. We refer to [16] for an overview of probability-constrained problems and list below some of our motivating settings and references:

Microgrid management: An electric power microgrid is a collection of intermittent renewable generator units, a conventional dispatchable diesel generator (or grid interconnection), and a battery energy storage system. The microgrid supplies electricity to a community in islanded mode, balancing fluctuating demand and supply. The operator achieves this by optimizing the use of the battery storage and the back-up dispatchable generator. Since perfect balancing is very expensive, it is common to allow for a small frequency of blackouts, i.e. occurrences where demand outstrips supply. Mixed-integer linear programming approaches to this problem through approximating with more conservative convex constraints appear in [23, 34]. 2. 2.

Hydropower optimization: control of a hydropower dam with probabilistic constraints was discussed in [2]. Within this setup, the controller observes random inflows from precipitation, as well as fluctuating electricity prices. His objective is to control the downstream outflow from the dam to maximize profit from power sales, while ensuring a minimum dam capacity with high probability. Other related works are [37, 3, 32]. 3. 3.

Motion planning: finding the minimum-cost path for a robot from one location to another while avoiding colliding with objects that obstruct its path. Stochasticity in the environment implies that the robot motion is only partially controlled. Robust optimization that guarantees obstacle avoidance might be infeasible, making probabilistic constraints a viable alternative. Dynamic programming methods for unmanned aerial vehicles were introduced in [33] and the probabilistic-constrained motion of a robot was solved in [18].

Contribution. We consider state-dependent probabilistic constraints that are expressed through an expectation constraint at each system state. Our setup involves a continuous-state, continuous-time model with a discrete-step control, where the constraints are imposed step-by-step. Therefore, at each step and at each state, the controller must estimate which controls are admissible and then optimize over the latter. While this setting is simpler than global probabilistic constraints, it is much harder than unconstrained control, because a secondary numerical procedure is needed as part of dynamic programming to repeatedly compute the admissibility sets. The canonical setup involves finite-horizon control of a stochastic process described through a stochastic differential equation of Itô type. The solution paradigm involves the Bellman or Dynamic Programming equation (DPE), which works with discretized time-steps, but with a smooth spatial variable. In this context, we develop algorithms to solve stochastic optimal control problems with probabilistic constraints using regression Monte Carlo (RMC). To make this highly nontrivial extension to RMC, we investigate tools from machine learning (including support vector machines (SVM), Gaussian process (GP) regression, parametric density estimation, logistic regression and quantile regression) to statistically estimate the admissible set as a function of the system state. Our algorithm handles the two parts of the problem—the constraint estimation and the approximation of the conditional expectation—in parallel and with significantly lower simulation budget compared to a naive implementation.

After benchmarking the proposed approaches on two practical case-studies from energy battery management, our main finding is to recommend logistic regression and GP-smoothed probability estimation as the best procedures. These methods are stable, relatively fast and allow for a variety of further adjustments and speed-ups. In contrast, despite theoretical appeal, quantile regression and SVM are not well-suited for this task. On a higher level, our main take-away is that DPE-based stochastic control with probabilistic constraints (SCPC) is well within reach of cutting-edge RMC methods. Thus, it is now computationally feasible to tackle such problems, opening the door for new SCPC models and applications.

Related Models for Probabilistic Constraints. There is an extensive literature on one-period optimization with chance constraints and on global multi-period probabilistic constraints. For the one-period formulations, the most popular approach is to transform the problem into linear or non-linear programs over a set of scenarios [11, 12, 28, 31]. In particular, Monte Carlo scenarios as employed below are very common, but the typical setup involves a single optimization problem, while we face an infinite family of them indexed by the system state $\mathbf{x}$ and the time-step $n$. Global probabilistic constraints in multi-period settings are tackled from multiple perspectives. The dynamic programming method [30] incorporates the constraint into the objective function via a Lagrange multiplier. The solution is then obtained by iteratively solving for the optimal control and the Lagrange multiplier. However, the solution is sub-optimal due to the duality gap. Mixed-integer linear programming [1] works by linearizing the constraints and requires to discretize the state space. An alternative [6, 7] is to transform into a static problem which is however computationally feasible only under strict assumptions on the system dynamics and noise distribution. Yet another option is the stochastic viability approach [2, 15] that focuses on maximizing the probability of being admissible, which is defined both in terms of profit targets and satisfying constraints at every time step.

Compared to above, our model of a multi-period optimization with one-step probabilistic constraints applied at each time step is new in the literature. To our knowledge, the closest setup is studied recently in [19] to compute the hedging price of a portfolio whose risk is defined in terms of its future value with respect to a set of stochastic benchmarks. Besides a local probabilistic constraint, the authors also provide dynamic programming equations for multi-period constraints. However, their solution is driven by very specific loss functions and state processes. In contrast, we develop general purpose numerical schemes using statistical learning methods.

2 Problem formulation

We study numerical resolution of stochastic control problems on finite horizon $[0,T]$ with local implicit constraints, specifically we work with constraints defined through probabilistic conditions on the controlled state. Let $(\mathbf{X}(t))_{t\geq 0}\in\mathcal{X}\subset\mathbb{R}^{d}$ be a continuous time controlled Markov process adapted to a given filtration $(\mathcal{F}_{t})$. The control is an $(\mathcal{F}_{t})$-adapted process $(u(t))_{t\geq 0}$, taking values in $\mathcal{W}\subset\mathbb{R}$. We further assume that control decisions are made at discrete epochs $\{t_{0},t_{1},\ldots,t_{N}=T\}$; between time-steps the value of $u(t)$ remains constant. Thus, the control process is piecewise-wise constant and càdlàg (right-continuous with left limits), and will be alternatively represented as $u(t)=\sum_{n=0}^{N-1}u_{n}\mathbbm{1}_{[t_{n},t_{n+1})}(t)$.

Our setup is essentially discrete-time as far as control is concerned, but we introduce the continuous-time system state because admissibility of actions depends on the trajectory of $\mathbf{X}(t)$ between control points. In our motivating examples, the dynamics of the system is described by a stochastic differential equation:

[TABLE]

where $(B(t))$ is a $m-$dimensional Brownian motion generating the filtration $(\mathcal{F}_{t})$ and $b:\mathbb{R}_{+}\times\mathcal{X}\times\mathcal{W}\rightarrow\mathbb{R}^{d}$ and $\sigma:\mathbb{R}_{+}\times\mathcal{X}\times\mathcal{W}\rightarrow\mathbb{R}^{d\times m}$ are measurable functions, such that a unique (weak) solution exists for admissible controls defined below and takes values in $\mathcal{X}$. For the convenience of notation, we will write $\mathbf{X}_{n}$ for $\mathbf{X}(t_{n})$, $n=0,\ldots,N$, and will not indicate explicitly the dependence on the control $(u(t))$ if the latter is clear from the context.

Admissibility is defined in feedback form via

[TABLE]

for deterministic functions $P_{k}:\mathcal{X}\times\mathcal{W}\to\mathbb{R}$ and given subsets $\mathcal{A}_{k}\subset\mathbb{R}$. For the rest of the article we assume $P_{n}(\mathbf{X}_{n},u_{n})$ and $\mathcal{A}_{n}$ to be of the form

[TABLE]

where $\mathcal{G}_{n}$ is a functional defined on trajectories of $(\mathbf{X}(s))$ over the time interval ${s\in[t_{n},t_{n+1})}$. In other words, we target the set of controls such that the conditional probability of the “failure” functional $\mathcal{G}_{n}(\cdot)$ of $\mathbf{X}$ being greater than zero is bounded by a threshold $p$, i.e.

[TABLE]

The threshold $p$ in equation (3) is interpreted as relaxing the strong constraint $\mathcal{G}_{n}(\cdot)\leq 0$ which may not be appropriate in a stochastic environment. Typical values of $p$ would generally be small ($p\approx 0.05$). We assume in the following that at least one admissible control exists at any state and at any time, hence $\mathcal{U}_{n:N}$ is determined by $\mathcal{U}_{k}(\mathbf{x})=\mathcal{U}_{k:k}(\mathbf{x})$, $k=n,\ldots,N-1$, the non-empty sets of admissible controls satisfying the constraints at decision epochs $t_{k}$ conditional on $\mathbf{X}_{k}=\mathbf{x}$.

The controller must optimize as well as evaluate the feasibility of proposed actions, with the performance criterion of the form

[TABLE]

where $W(\cdot)$ represents the terminal penalty and $\pi_{t}(\cdot,\cdot)$ the running cost. We re-write (4) in terms of the corresponding dynamic programming equation at step $n$:

[TABLE]

Above $\mathcal{C}_{n}(\mathbf{X}_{n},u)$ is the continuation value, i.e. reward-to-go plus expectation of future rewards, from using the control $u$ over $[t_{n},t_{n+1})$. Moreover, given the state $\mathbf{X}_{n}$, we say that $u^{*}\in\mathcal{U}_{n}(\mathbf{X}_{n})$ is an optimal control if $V_{n}(\mathbf{X}_{n})=\mathcal{C}_{n}(\mathbf{X}_{n},u^{*})$. Since the admissible set $\mathcal{U}_{n}(\mathbf{X}_{n})$ is both time and state dependent, we need to estimate the continuation value $\mathcal{C}_{n}(\cdot,\cdot)$ and the admissible control set $\mathcal{U}_{n}(\cdot)$ at every time step. This is the major distinction from the standard scenario approach [28] in chance-constrained optimization where there is only a single problem to optimize over a fixed $\mathcal{U}$, but no further indexing by $x$ and by $n$. The latter require a combination interpolation and optimization as part of the solution.

Alternative Formulation of Admissibility. We denote by $G_{n}(\mathbf{X}_{n},u_{n})$ as the regular conditional distribution [20] of the functional $\mathcal{G}_{n}(\cdot)$ given $(\mathbf{X}_{n},u_{n})$:

[TABLE]

where $\mathcal{L}(\cdot|X_{n},u_{n})$ stands for a conditional law. When writing \mathbb{P}\big{(}G_{n}(\mathbf{X}_{n},u_{n})>z\big{)} or \mathbb{E}\big{[}g\big{(}G_{n}(\mathbf{X}_{n},u_{n})\big{)}\big{]} we mean the probability or the expectation with respect to this conditional distribution.

We may rewrite equation (2) through the corresponding $(1-p)^{th}$ quantile $q(\mathbf{X}_{n},u_{n})$ of $G_{n}(\mathbf{X}_{n},u_{n})$:

[TABLE]

Then using

[TABLE]

we can set $P^{\prime}_{n}:=q_{n}$ and $\tilde{\mathcal{A}}=(-\infty,0]$ in (1). We will exploit this equivalence to propose quantile-based methods (Section 4) for the admissible set.

Remark 2.1**.**

Assuming a one dimensional control $u_{n}\in\mathcal{W}\subset\mathbb{R}$, and the probability $p_{n}(\mathbf{X}_{n},u_{n})$ monotonically decreasing in $u_{n}$, estimating the admissible set $\mathcal{U}_{n}(\mathbf{X}_{n})$ is equivalent to estimating the minimum admissible control

[TABLE]

*The corresponding admissible set is $\mathcal{U}_{n}(\mathbf{X}_{n})=\{u\in\mathcal{W}:u\geq u^{\min}_{n}(\mathbf{X}_{n})\}$. *

Remark 2.2**.**

A more general version are implicit constraints of the form

[TABLE]

*for a function $g:\mathbb{R}\to\mathbb{R}$. Even more abstractly, we can think of a generic implicit map $P_{k}(\cdot,\cdot)$ in (1) that defines $\mathcal{U}_{n:N}(\mathcal{X}_{n})$, with the idea that inverting this map is numerically nontrivial (i.e. $P_{k}$ is expensive to evaluate) and hence a priori it is not clear which controls satisfy constraints and which do not. *

Remark 2.3**.**

Equation (8) describes admissible controls $u$ for a given state $\mathbf{x}$. The “dual” perspective is to consider the set of states $\mathcal{X}^{a}_{n}(u)\subset\mathcal{X}$ for which a given control $u$ is admissible:

[TABLE]

*Often the cardinality of $\mathcal{X}$ is infinite, while the control space $\mathcal{W}$ is finite, so that enumerating (9) over $u\in\mathcal{W}$ is considerably easier than enumerating the uncountable family of sets $\mathbf{x}\mapsto\mathcal{U}_{n}(\mathbf{x})$ in equation (3). Furthermore, if $u\mapsto p_{n}(\mathbf{x},u)$ is decreasing for all $\mathbf{x}\in\mathcal{X}$, then we obtain an ordering $\mathcal{X}^{a}_{n}(u_{1})\subseteq\mathcal{X}^{a}_{n}(u_{2})$ for $u_{1}\leq u_{2}$. The latter nesting feature corresponds to ranking the controls in terms of their “riskiness” with respect to $G_{n}$, so that the safest control will have a very large $\mathcal{X}^{a}_{n}(u)$ (possibly all of $\mathcal{X}$), while the riskiest control will have a very small admissibility domain. *

2.1 Regression Monte Carlo

In this article we focus on simulation-based techniques to solve (4). The overall framework is based on solving equation (5) through backward induction on $n=N-1,N-2,\ldots$, replacing the true $V_{n}(\mathbf{x})$ with an estimate $\hat{V}_{n}(\mathbf{x})$. Since neither the conditional expectation, nor the admissibility constraint are generally available explicitly, those terms must also be replaced with their estimated counterparts. As a result, we work with the approximate Dynamic Programming recursion

[TABLE]

Above, $\hat{\mathbb{E}}$ is the approximate projection operator and the set of admissible controls $\hat{\mathcal{U}}_{n}$ is also approximated via either $\hat{p}_{n}(\cdot,\cdot)$, i.e., \hat{\mathcal{U}}_{n}(\mathbf{X}_{n}):=\big{\{}u:\hat{p}_{n}(\mathbf{X}_{n},u)<p\big{\}}, or $\hat{q}_{n}(\cdot,\cdot)$, i.e., \hat{\mathcal{U}}_{n}(\mathbf{X}_{n})=\big{\{}u:\hat{q}_{n}(\mathbf{X}_{n},u)\leq 0\big{\}}, see (8). The estimated optimal control $\hat{u}_{n}\in\hat{\mathcal{U}}_{n}(\mathbf{X}_{n})$ satisfies $\hat{V}_{n}(\mathbf{X}_{n})=\hat{\mathcal{C}}_{n}(\mathbf{X}_{n},\hat{u}_{n}).$

The key idea underlying our algorithm and defining the Regression Monte Carlo paradigm is that $\hat{\mathbb{E}}$ and $\hat{\mathcal{U}}$ are implemented through empirical regressions based on Monte Carlo simulations. In other words, we construct random, probabilistically defined approximations based on realized paths of $\mathbf{X}$. This philosophy allows to simultaneously handle the numerical integration (against the stochastic shocks in $\mathbf{X}$) and the numerical interpolation (defining $\hat{V}_{n}(x)$ for arbitrary $x$) necessary to solve (10).

To understand RMC, recall that specifying $\hat{\mathbb{E}}$ is equivalent to approximating the conditional expectation map (\mathbf{x},u)\mapsto\mathbb{E}[\psi\big{(}(\mathbf{X}(s))_{s\in[t_{n},t_{n+1}]}\big{)}|\mathbf{X}_{n}=\mathbf{x},u_{n}=u]=:f(\mathbf{x},u) where we specifically substitute

[TABLE]

To do so, we consider a dataset consisting of inputs $(\mathbf{x}^{1}_{n},u^{1}_{n}),\ldots,(\mathbf{x}^{M_{c}}_{n},u^{M_{c}}_{n})$ and the corresponding pathwise realizations $y^{1},\ldots,y^{M_{c}}$ with y^{j}=\psi\big{(}(\mathbf{x}(s))^{j}_{s\in[t_{n},t_{n+1}]}\big{)}, where $(\mathbf{x}(s))^{j}_{s\in[t_{n},t_{n+1}]}$ is an independent draw from the distribution of the process $(\mathbf{X}(s))_{s\in[t_{n},t_{n+1}]}|(\mathbf{x}^{j}_{n},u^{j}_{n})$. Then we use the training set $\{\mathbf{x}^{j}_{n},u^{j}_{n},y^{j}\}_{j=1}^{M_{c}}$ to compute $\hat{f}$, an estimator of $f$, via regression.

Similarly, estimating $\mathcal{U}_{n}$ is equivalent to learning the conditional probability map $p_{n}(\mathbf{x},u)$ (or the conditional quantile map ${q}_{n}(\mathbf{x},u)$ in (3)) and then comparing to the threshold value $p$ (zero, respectively). This statistical task, whose marriage with RMC is our central contribution, is discussed in Section 4.

The technique of using regressions for the approximation of the continuation value was developed in the celebrated works by [24] and [36] in the context of American option pricing and further enhanced in [10, 25]. This was extended for storage problems and controlled state process in [13, 8, 9, 4, 26]. Among the approaches for approximating $f$ we mention [39] and [26] who exploit the structure of the problem to reduce the dimensionality of the regressions, [4, 5] who harness the distribution of process to reduce the variance of $\hat{f}$ and [39, 21, 26] who use non-parametric regression methods for $\hat{f}$. Regression based approach has also been discussed in the context of stochastic dual dynamic programming for solving high dimensional storage problems in [38].

In contrast to the above well-developed literature, very little exists about estimating the set of admissible controls $\hat{\mathcal{U}}_{n}(\mathbf{X}_{n})$, which requires approximating $p(\mathbf{x},u)$ (or $q(\mathbf{x},u)$) in equation (3). There are results about learning a single global admissibility set $\hat{\mathcal{U}}$ in a one-period setup, but those approaches do not transfer to our context of state- and time-dependent admissibility constraints. A naive approach is to estimate $\hat{\mathcal{U}}_{n}(\mathbf{X}_{n})$ for every state realized during the backward induction through nested Monte Carlo. Namely for each pair $(\mathbf{x},u)$ encountered, we may estimate the probability of violating the constraint by simulating $M_{b}$ samples from the conditional distribution $G_{n}(\mathbf{x},u)$ as $\{g_{n}^{b}(\mathbf{x},u)\}_{b=1}^{M_{b}}$. We then set $u\in\hat{\mathcal{U}}_{n}(\mathbf{x})$ if $\bar{p}_{n}(\mathbf{x},u)<p$, where

[TABLE]

is the empirical probability. Although simple to implement, this Nested Monte Carlo (NMC) method is computationally intractable for even the easiest problems. As an example, a typical RMC scheme employs $M_{c}\approx 100,000$ and assuming $M_{b}=1000$ for inner simulations, which is necessary for good estimates of small probabilities $p\leq 0.1$, would require $10^{8}$ simulation budget at every time-step to implement NMC. Note furthermore that NMC returns only the local estimates $\bar{p}(\mathbf{x},u)$; no functional estimate of $\mathcal{U}_{n}(\mathbf{x})$ or $\mathcal{X}^{a}_{n}(u)$ is provided for an arbitrary $\mathbf{x}$ or $u$, respectively. As a result, any out-of-sample evaluation (i.e. on a future sample path of $\mathbf{X}$) requires further inner simulations, making this implementation even more computationally prohibitive.

An important challenge in using $\hat{\mathcal{U}}$ is verifying admissibility. Since we are employing Monte Carlo samples to decide whether $u$ is admissible at $\mathbf{x}$, this is a probabilistic statement and admissibility can never be guaranteed 100%. We may use statistical theories to quantify the accuracy of estimators of $\mathcal{U}$, for example, by applying Central Limit Theorem tools for the estimator $\bar{p}(\mathbf{x},u)$ of the true $p(\mathbf{x},u)$. In particular, we develop tools based on confidence intervals in order to make statements (with asymptotic guarantee) such as “$u\in\mathcal{U}$ with 95% confidence” (equivalent to $p(\mathbf{x},u)<p$ with 95% probability conditional on the data collected). Achieving reasonable confidence levels calls for “conservative” estimators of $\hat{\mathcal{U}}$. As we show, not doing so can make learning $\mathcal{U}$ highly unreliable, frequently causing decisions that are inadmissible with respect the imposed probability constraint. Thus, the related construction of $\hat{\mathcal{U}}^{(\rho)}$ with specified confidence level $\rho$ is a running theme in Section 4.

2.2 Motivation: controlling blackout probability in a microgrid

To make our presentation concrete, we illustrate the framework of (5) by formalizing the motivating application from microgrid management. A microgrid comprises renewable and traditional generation sources, along with a medium of storage, designed and managed to provide electrical power to a community in a decentralized way. We consider a system composed of a dispatchable diesel generator, a renewable energy source and an electric battery storage. The microgrid topology is illustrated in the left panel of Figure 1 and is same as the example discussed in [26].

In this context, the state variables are $\mathbf{X}(t)=(L(t),I(t),C(t))$, where $L(t)$ is the net demand (demand net of renewable generation), $I(t)\in[0,I_{\max}]$ is the state of charge of the battery, referred to as “the inventory”, and $C(t)\in\{0,1\}$ is the state of the diesel generator. $C(t)=0$ refers to diesel being OFF and $C(t)=1$ implies ON. The controller is in charge of the diesel through the control $u(t)$, which indicates the power output of the unit. We assume, for clarity of exposition, that the net demand $L(t)$ is an exogenous process, while $I(t)$ is controlled. We reiterate that the control decisions are made at discrete epochs $\{t_{0},t_{1},\ldots,t_{N-1}\}$, however these decisions affect the state of the system continuously. The choice of $u(t_{n})\equiv u_{n}$ at time $t_{n}$ is based on minimizing the cost of running the microgrid, as well as controlling the probability of a blackout (i.e. failing to match the net demand) during $[t_{n},t_{n+1})$. The blackout is described through the imbalance process $S(s):=L(s)-u_{n}-B(s)$, $\forall s\in[t_{n},t_{n+1})$, representing the difference between the demand and supply, while the diesel output is held constant over the time step. The power output from the battery is a deterministic function of net demand, inventory and the control, $B(s)=\varphi(L(s),I(s),u_{n})$ constrained by the physical limitations of the battery. $B(s)>0$ implies supply of power from the battery and $B(s)<0$ implies battery charging. The set of admissible controls is thus:

[TABLE]

Thus in the context of microgrid, the conditional distribution $G_{n}$ of equation (6) and the corresponding $p_{n}(L_{n},I_{n},C_{n})$ are:

[TABLE]

Because $p_{n}$ is not (in general) available analytically, the admissibility condition $p_{n}(L_{n},I_{n},C_{n})<p$ is implicit. Recall that we denote by $\mathcal{W}=0\cup[\underline{u},\bar{u}]$ the unconstrained control set. We assume that $u(t)=0$ means that the diesel is OFF, while $u(t)>0$ means that it is ON, and at output level $u(t)$. Thus, we define $C(s)=\mathbbm{1}_{\{u_{n}>0\}}$ $\forall s\in(t_{n},t_{n+1}]$ with the time interval left-open in order to allow for identification of switching on and off of the diesel generator at times $t_{n}$. Notice also that the process $C(t)$ does not satisfy the controlled diffusive dynamics, but this slight extension of the framework does not impact on the methods and results presented. We then look at the following formulation of the general problem:

[TABLE]

where $\Delta t_{k}=t_{k+1}-t_{k}$, $\rho(u_{k})$ is the instantaneous cost of running the diesel generator with power output $u_{k}$ and $\mathcal{K}$ is the cost of switching it ON. We assume zero cost to turn the generator off. The DPE corresponding to (14) is the same as in (5) with the integral running cost $\int_{t_{n}}^{t_{n+1}}\pi_{s}(\mathbf{X}(s),u_{n})\mathrm{d}s$ replaced by

[TABLE]

Remark 2.4**.**

*The admissible set $\mathcal{U}\subseteq\mathcal{W}$ for this problem has the special structure: if $u\in\mathcal{U}(\mathbf{x})$, then $\forall\ \mathcal{W}\ni\tilde{u}>u,\tilde{u}\in\mathcal{U}(\mathbf{x})$. Hence, we may represent $\mathcal{U}(\mathbf{x})=[u^{\min}_{n}(\mathbf{x}),\bar{u}]\cap\mathcal{W}$ in terms of the minimal admissible diesel output $u^{\min}_{n}(\mathbf{x})$. Conversely, the admissibility domains for a fixed $u\in\mathcal{W}$ are nested: if $u_{1}\leq u_{2}$ then $\mathcal{X}^{a}_{n}(u_{1})\subseteq\mathcal{X}^{a}_{n}(u_{2})$. This suggests to compute $\mathcal{X}^{a}_{n}(u)$ sequentially as $u$ is increased and then invert to get $\mathcal{U}(\mathbf{x})$. *

To visualize the minimum admissible control $u^{\min}_{n}(\mathbf{x})$, the right panel of Figure 1 presents the map $\mathbf{x}\rightarrow u^{\min}_{n}(\mathbf{x})$ under a constraint of $p=0.01$ probability of blackout. We also present a path for $(L(t),I(t),C(t))_{t\geq 0}$ using a myopic strategy where the controller employs the minimum admissible control at each point, $u_{n}:={u}^{\min}_{n}(L_{n},I_{n},C_{n})\;\forall n$. Notice how for the most part, $u^{\min}_{n}(\cdot)=0$ so that $\mathcal{U}_{n}(\cdot)=\mathcal{W}$ and the blackout constraint is not binding. This is not surprising, as blackouts are only possible when $L(t)\gg 0$ is strongly positive and the battery is close to empty, $I(t)\simeq 0$. Thus, except for the lower-right corner, any control is admissible. As a result, only a small subset of the domain $\mathcal{X}$ actually requires additional effort to estimate the admissible set $\mathcal{U}(\mathbf{x})$. In our experience this structure, where the constraint is not necessarily binding and where we mostly perform unconstrained optimization, is quite common.

3 Dynamic emulation algorithm

In this section we present our Dynamic emulation algorithm which provides approximation for the admissible set $\hat{\mathcal{U}}_{n}(\cdot)$ and the continuation value function $\hat{\mathcal{C}}_{n}(\cdot,\cdot)$. The crux of the algorithm are the following two steps, implemented in parallel at every time-step:

[TABLE]

To estimate $\hat{\mathcal{C}}_{n}(\cdot,\cdot)$’s and $\hat{\mathcal{U}}_{n}(\cdot)$’s, we proceed iteratively backward in time starting with known terminal condition $W(\mathbf{X})$ and sequentially estimate $\hat{\mathcal{U}}_{n}$ and $\hat{\mathcal{C}}_{n}$ for $n=N-1,\ldots,0$. Assuming we have estimated $\hat{\mathcal{U}}_{n+1},\ldots,\hat{\mathcal{U}}_{N-1}$ and $\hat{\mathcal{C}}_{n+1},\ldots,\hat{\mathcal{C}}_{N-1}$, we first explain the estimation procedure for $\hat{\mathcal{U}}_{n}$ and $\hat{\mathcal{C}}_{n}$. This corresponds to a fit task. In the subsequent backward recursion at step $n-1$ we also need the predict task to actually evaluate $\hat{V}_{n}(\mathbf{X}_{n})$ which requires evaluating $\hat{\mathcal{C}}_{n}(\cdot)$ at new (“out-of-sample”) inputs $\mathbf{X}_{n},u_{n}$ which of course do not coincide with the training inputs $(\mathbf{x}^{1}_{n},u^{1}_{n}),\,\dots,\,(\mathbf{x}^{M_{c}}_{n},u^{M_{c}}_{n})$.

3.1 Estimating the set of admissible controls

To estimate the set of admissible controls $\hat{\mathcal{U}}_{n}(\cdot)$ at time-step $n$, we choose design $\mathcal{D}^{a}_{n}:=(\mathbf{x}_{n}^{i},u_{n}^{i},i=1,\ldots,M_{a})$ and simulate trajectories of the state process $(\mathbf{X}(s))^{i}_{s\in[t_{n},t_{n+1})}$ starting from $\mathbf{X}^{i}(t_{n})=\mathbf{x}_{n}^{i}$ and driven by control $u_{n}^{i}$. To evaluate the functional \mathcal{G}\big{(}(\mathbf{X}(s))^{i}_{s\in[t_{n},t_{n+1})}\big{)}, we discretize the time interval $[t_{n},t_{n+1})$ into $K$ finer sub-steps with $\Delta n_{k}:=t_{n_{(k+1)}}-t_{n_{k}}$ and define the discrete trajectory $\mathbf{x}_{n}^{i}=\mathbf{x}_{n_{0}}^{i},\mathbf{x}_{n_{1}}^{i},\ldots,\mathbf{x}_{n_{(K-1)}}^{i},\mathbf{x}_{n_{K}}^{i}$. We then record

[TABLE]

where, formally, we extend $(\mathbf{x}_{n_{k}}^{i})_{k\in\{0,\ldots,K-1\}}$ to a piecewise constant trajectory on $[t_{n},t_{n+1})$.

Analogous to standard RMC, we now select an approximation space $\mathcal{H}^{a}_{n}$ to estimate the probability $\hat{p}_{n}$ or the quantile $\hat{q}_{n}$, using the loss function $\mathcal{L}^{a}_{n}$ and apply empirical projection:

[TABLE]

See Section 4 for concrete examples of $\mathcal{H}^{a}$ and $\mathcal{L}^{a}$. Note that the approximations $\hat{p}_{n}$ and $\hat{q}_{n}$ must be trained on joint state-control datasets $\{\mathbf{x}_{n}^{i},u_{n}^{i},w^{i}_{n}\}_{i=1}^{M_{a}}$ with $w^{i}_{n}$ dependent on the method of choice and moreover yield random estimators ($\hat{p}_{n}$ is a random variable).

Using the distribution of $\hat{p}_{n}(\mathbf{x},u)$ we may obtain a more conservative estimator that provides better guarantees on the ultimate admissibility of $(\mathbf{x},u)$. As a motivation, recall the NMC estimator $\bar{p}_{n}(\mathbf{x},u)$ from (11); for reasonably large $M_{b}\gg 20$, the distribution of $\bar{p}_{n}(\mathbf{x},u)$ is approximately Gaussian with mean $p_{n}(\mathbf{x},u)$ and variance $\frac{p_{n}(\mathbf{x},u)(1-p_{n}(\mathbf{x},u))}{M_{b}}$. Defining

[TABLE]

where $z_{\rho}$ is the standard normal quantile at level $\rho$ and $\xi^{(\rho)}_{n}(\mathbf{x},u)$ represents a “safe” margin of error for $\bar{p}_{n}$ at confidence level $\rho$. The corresponding approximate admissible set with confidence $\rho$ is

[TABLE]

More generally, we set the admissible set for a site $\mathbf{x}\in\mathcal{X}$ to

[TABLE]

where $\xi_{n}(\mathbf{x},u)$ ensures “stronger” guarantee for the admissibility of $u$ at $\mathbf{x}$. The margin of estimation error can also be fixed, $\xi_{n}(\mathbf{x},u)=c\;\forall(\mathbf{x},u)\in\mathcal{X}\times\mathcal{W}$, which can be applied when the sampling distribution of $\hat{p}_{n}(\mathbf{x},u)$ is unknown. The corresponding admissible set

[TABLE]

is equivalent to estimating $\hat{\mathcal{U}}_{n}^{\xi=0}(\mathbf{x})$ with a shifted lower probability threshold $p-c$. To simplify notation, we use $\hat{\mathcal{U}}_{n}(\mathbf{x})$ to denote the unadjusted admissible set, $\hat{\mathcal{U}}_{n}(\mathbf{x}):=\hat{\mathcal{U}}_{n}^{\xi=0}(\mathbf{x})$ in the context of NMC. Analogously, we can adjust equations (16)-(17) based on learning the quantile $q_{n}(\mathbf{x},u)$ to add a margin of error, $\hat{\mathcal{U}}_{n}^{\xi}(\mathbf{x})=\{u:\hat{q}_{n}(\mathbf{x},u)+\xi_{n}(\mathbf{x},u)\leq 0\}$.

Remark 3.1**.**

*In the context of one-step static optimization, conservative estimates for the admissible set were explored in [12]. The idea is to consider $u$ to be admissible iff $\overline{p}(\mathbf{x},u)=0$, i.e. there are no failures observed in all of the simulated samples. The choice of the sample size $M_{b}$ is then determined using Chebyshev inequality that provides a theoretical guarantee on $\mathbb{P}(p(\mathbf{x},u)>0\,|\,\overline{p}(\mathbf{x},u;M_{b})=0)$. This approach is similar in our setup to a violation margin $c=p$ in (22). However, the guarantee on $p(\mathbf{x},u)$ is only applicable locally at the design sites. We are not aware of any tools that would allow non-parametric guarantee for any $\mathbf{x}\in\mathcal{X}$. The regression based approach offers a model-based guarantee for admissibility of $u$ at any $\mathbf{x}$ by setting the parameter $\rho$ in (20) or $c$ in (22). *

3.2 Estimating the continuation value

To estimate the continuation value $\mathcal{C}_{n}(\cdot,\cdot)$, we choose a simulation design $\mathcal{D}^{c}_{n}:=(\mathbf{x}_{n}^{j},u_{n}^{j},j=1\ldots,M_{c})$ (which could be independent or equivalent to $\mathcal{D}^{a}_{n}$) and generate one-step paths for the state process $(\mathbf{X}(s))^{j}_{s\in[t_{n},t_{n+1})}$ starting from $\mathbf{X}^{j}(t_{n})=\mathbf{x}_{n}^{j}$ and driven by control $u_{n}^{j}$, comprising again finer sub-steps $\mathbf{x}_{n}^{j}=\mathbf{x}_{n_{0}}^{j},\mathbf{x}_{n_{1}}^{j},\ldots,\mathbf{x}_{n_{(K-1)}}^{j},\mathbf{x}_{n_{K}}^{j}$ (in principle the sub-steps could differ from the time discretization for $\hat{\mathcal{U}}_{n}$). Next, we compute the pathwise cost $y_{n}^{j}$:

[TABLE]

and we replace the time integral in (5) with a discrete sum over $t_{n_{k}}$’s. At the key step, we project $\{y_{n}^{j}\}_{j=1}^{M_{c}}$ onto an approximation space $\mathcal{H}^{c}_{n}$ to evaluate the continuation value $\mathcal{C}_{n}(\cdot,\cdot)$:

[TABLE]

The design sites $\{\mathbf{x}_{n}^{j},u_{n}^{j}\}_{j=1}^{M_{c}}$ could be same or different from those used for learning the admissible sets in the previous subsection. Two standard approximation spaces $\mathcal{H}^{c}_{n}$ used in this context are: global polynomial approximation and piecewise continuous approximation.

Remark 3.2**.**

*In the microgrid example of Section 2.2 the running cost over $[n,n+1)$ is known once the control $u_{n}$ is chosen. Thus it can be taken outside the conditional expectation and the data to be regressed is simply $y^{j}=v_{n+1}^{j}$. *

Global polynomial approximation:

This is a classical regression framework with polynomial bases $\phi_{k}(\cdot,\cdot)$ and $\hat{\mathcal{C}}_{n}^{\bm{\alpha}}(\mathbf{x},u):=\sum_{k}\alpha_{k}\phi_{k}(\mathbf{x},u)$. The coefficients $\bm{\alpha}$ are fitted via

[TABLE]

As an illustration, for the microgrid example of Section 2.2 we construct a quadratic polynomial approximation when diesel generator is ON, $u>0$, using $10$ bases $\{1,L,I,u,L^{2},I^{2},u^{2},LI,Iu,LI\}$ and a separate quadratic approximation with the $6$ basis functions $\{1,L,I,L^{2},I^{2},LI,LI\}$ when diesel generator is OFF, $u=0$. Polynomial approximation is easy to implement but typically requires many degrees of freedom (lots of $\phi$’s) to properly capture the shape of $\mathcal{C}$ and can be empirically unstable, especially if there are sharp changes in the underlying function (see for example [26, 21]).

Piecewise continuous approximation:

This is a state-of-art tool in low dimensions, $d\leq 3$. The main idea is to employ polynomial regression in a single dimension and extend to the other dimensions via linear interpolation. As an example, for the microgrid with diesel generator ON, we have three dimensions $(L,I,u)$. We discretize inventory $I$ as $\{I^{0},I^{1},\ldots,I^{M_{I}}\}$ and control $u$ as $\{u^{1},u^{2},\ldots,u^{M_{u}}\}$ and fit independent cubic polynomials in $L$ for each pair $(I^{l},u^{e})$ with $l\in\{0,1,\ldots,M_{I}\}$ and $e\in\{1,\ldots,M_{u}\}$, i.e., $f_{n}^{l,e}(L)=\sum_{k}\alpha_{k}^{l,e}\phi_{k}(L)$. For any $(I,u)\in[I^{l},I^{l+1}]\times[u^{e},u^{e+1}]$ we then provide the interpolated approximation $\hat{\mathcal{C}}_{n}(L,I,u)$ as

[TABLE]

Nonparametric approximation:

Further alternatives for $\mathcal{H}_{n}^{c}$ can be found in [26] who used Gaussian process regression and [21, 22] who used local polynomial regression. For semi-parametric approximation, [10] developed piecewise multivariate linear regression.

There are several possibilities for choosing the designs $\mathcal{D}^{a}_{\cdot}$ and $\mathcal{D}^{c}_{\cdot}$, see [26] for a detailed discussion of different regression designs $\mathcal{D}^{c}_{\cdot}$ and their impact on the quality of the final solution.

3.3 Evaluation

We analyze the quality of the solution by computing three quantities on the out-of-sample dataset:

• •

estimate of the value function $V_{0}(\mathbf{x}_{0})$ at $t=0$ and state $\mathbf{x}_{0}$;

• •

empirical frequency of inadmissible decisions on the controlled trajectories $\mathbf{x}^{\hat{u}}_{\cdot}$;

• •

statistical test for the realized number of constraint violations (blackouts for the microgrid).

Good solutions should minimize costs and not apply inadmissible controls. However, since we employ empirical estimators, $\mathcal{U}$ is never known with certainty and we must handle the possibility that constraints are violated with probability more than $p$. In turn this leads to the trade-off between complying with (1) and optimizing costs. Similar treatment of constraints in the context of sample average approximation of probabilistic constrained optimization problems have been discussed in [29, 27]. Moreover, our framework implies that the whole algorithm is stochastic: multiple runs will lead to different results since both $\hat{p}_{n}$ and $\hat{\mathcal{C}}_{n}$ are impacted by the random samples $y^{j}_{n}$ and $w_{n}^{i}$.

Estimate of the value function: We evaluate the value function $\hat{V}_{0}(\mathbf{x}_{0})$ at time $t_{0}=0$ and state $\mathbf{x}_{0}$ using $M^{\prime}$ out-of-sample paths $(\mathbf{x}_{0:N}^{\hat{u},m^{\prime}}),m^{\prime}=1,\ldots,M^{\prime}$. Each trajectory $(\mathbf{x}_{0:N}^{\hat{u},m^{\prime}})$ is generated by applying the estimated optimal control $\hat{u}_{0:N-1}$ based on the continuation value and admissible sets $(\hat{\mathcal{C}}_{n},\hat{\mathcal{U}}_{n})_{n=0}^{N-1}$ leading to the realized pathwise cost

[TABLE]

The resulting empirical Monte Carlo estimate is

[TABLE]

and represents an unbiased estimation of the value of the control policy and an asymptotic upper bound estimation of the value function, provided all controls used are admissible.

Empirical frequency of inadmissible decisions on the controlled trajectories: For the $M^{\prime}$ out-of-sample paths, we compare the estimated optimal control $\{\hat{u}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})\}_{n=m=1}^{N-1,M^{\prime}}$ against the minimum admissible control $\{u^{\min}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})\}_{n=m=1}^{N-1,M^{\prime}}$ assumed for a second to be known. Namely, for each path we compute the number of inadmissible decisions $w_{0}(\mathbf{x}_{0:N}^{\hat{u},m^{\prime}})$ and the respective empirical frequency $w_{freq}$ as:

[TABLE]

respectively. We employ these metrics in Section 5, where a “gold standard” $\{u^{\min}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})\}_{n=m=1}^{N-1,M^{\prime}}$ is obtained by brute force, utilizing a simulation budget $10^{5}$ larger than for the actual methods we are comparing. Empirical gold standard is a common technique when analytical benchmark is unavailable, see e.g. [14]. A good estimation method should yield $w_{freq}\simeq 0$. However, if $w_{freq}>0$, a controller can choose margin of error $\xi_{n}$ (Equation (21)) to reduce $w_{freq}$ at the expense of higher cost $v_{0}(\mathbf{x}_{0:N}^{\hat{u},m^{\prime}})$.

Remark 3.3**.**

*In Equation (28) we exploit the structure of the admissible set for the microgrid control (cf. Remark 2.4). Generally, $w_{0}(\mathbf{x}_{0:N}^{\hat{u},m^{\prime}}):=\sum_{n}\mathbbm{1}_{\hat{u}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})\notin\mathcal{U}_{n}(\mathbf{x}_{n})}$, where $\mathcal{U}_{n}(\mathbf{x}_{n})$ is either known in closed form or computed via empirical gold standard. The general setting, without any assumptions on the structure of $\mathcal{U}_{n}(\mathbf{x}_{n})$, is computationally very expensive and beyond the scope of this work. *

Statistical test: Next we propose statistical tests using the controlled trajectories to validate different methods for admissible set estimation. Such a test is essential to affirm the use of a numerical scheme for $\mathcal{U}_{n}$ in the absence of a benchmark. As an example, in the context of microgrid we want to test the null hypothesis $H_{0}$ that the realized probability of blackouts is bounded to the required level against the alternative $H_{1}$ that their probability is too high. Let

[TABLE]

Ignoring the correlation due to the temporal dependence in $\mathbf{x}_{n}$, we assume that $B_{n}^{m^{\prime}}\sim\text{Bernoulli}(\tilde{p})$ are i.i.d. and $\tilde{p}$ represents the true (unknown) probability of blackout. We want to test:

[TABLE]

A common approach to such composite null hypothesis is to replace $H_{0}$ with a more conservative hypothesis $\tilde{p}=p$ leading to the test statistic

[TABLE]

Hence, $H_{0}$ is rejected at a confidence level $\alpha$ if $\mathcal{T}>z_{\alpha}$ with $z_{\alpha}=\Phi^{-1}(\alpha)$, e.g. $z_{\alpha}=1.65$ for $\alpha=95\%$.

Remark 3.4**.**

*The above test assumes independence and identical distribution of $B_{n}^{m^{\prime}}$’s. In the context of the microgrid example, neither of the two assumptions are valid; $B_{n}^{m^{\prime}}$ has a different distribution because the state of the system affects the probability of a blackout, thus $\tilde{p}$ varies with $n,m^{\prime}$. Furthermore, $B_{n}^{m^{\prime}}$ are not independent as they are derived from a single, sequentially controlled trajectory. *

Remark 3.5**.**

In the microgrid setup, the blackout constraint is frequently not binding (the net demand is negative half of the time). Therefore, $\mathcal{T}$ as defined in equation (31) is most likely negative leading to accept the $H_{0}$ even when the method fails to choose the admissible control when the constraint is binding. We fix this by evaluating the sum only when the constraint is binding, i.e.

[TABLE]

*where $\hat{u}^{\min}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})$ represents model estimate of the true minimum admissible control $u^{\min}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})$. *

To wrap up this section, Algorithm 1 (dubbed Dynamic Emulation due to similarities with a related algorithm for unconstrained stochastic control from [26]) summarizes the overall sequence of steps. Lines 1-2 contain the parameters to the algorithm, Lines 3-8 (and 14-19) yield the stochastic simulator which generates designs and corresponding one-step paths. Line 10 (and again Line 20) computes pathwise one-step costs. Line 12 is the admissible set estimation. Line 13 is the estimation of the continuation value.

Algorithm 1 carries several advantages. First and foremost it is very general, and does not make any restrictions on the distribution $G_{n}(\mathbf{X}_{n},u)$ defining $\mathcal{U}_{n}$ or the form of the payoffs $\pi(\mathbf{x},u)$. Hence it can be generically applied across a wide spectrum of SCPC problems. Second, the same template (in particular based on having two independent sub-modules) accommodates a slew of techniques for learning $\mathcal{C}$ and $\mathcal{U}$ bringing plug-and-play functionality, such as straightforward switching from probability to quantile estimation. Third, it allows for computational savings through parallelizing the estimation of $\mathcal{U}$ and $\mathcal{C}$, or by re-using the same design and simulations $\mathcal{D}^{a}_{n}\equiv\mathcal{D}^{c}_{n}$ for the computation of the two sub-modules.

Remark 3.6**.**

The challenge of RMC methods is that the errors recursively propagate backward. As a result, poor estimation at one step can affect the overall quality of the solution. In our algorithm, the errors at every step occur due to:

• •

Approximation architecture $\mathcal{H}^{a}_{n}$ for $\hat{\mathcal{U}}_{n}$ $\Rightarrow$ Projection error in admissible control set estimation;

• •

Approximation architecture $\mathcal{H}^{c}_{n}$ for $\hat{\mathcal{C}}_{n}$ $\Rightarrow$ Projection error in estimating continuation value;

• •

Designs $\mathcal{D}^{a}_{n}$ and $\mathcal{D}^{c}_{n}\Rightarrow$ Finite-sample Monte Carlo errors (difference between empirical estimates and theoretical projection-based ones)

• •

Discretization of the time interval $[t_{n},t_{n+1})$ using $\Delta n_{k}$ $\Rightarrow$ Integration error in approximating the integral $\int_{t_{n}}^{t_{n+1}}\pi_{s}(\mathbf{X}(s),u)\mathrm{d}s$ and the admissible set $\mathcal{U}_{n}$.

• •

Numerical approximation of the solution of the controlled dynamics of $\mathbf{X}(t)$.

• •

Optimization errors in maximizing for $\hat{u}$ over $\hat{\mathcal{U}}$, especially when the control set $\mathcal{W}$ is continuous.

4 Admissible set estimation

In this section we propose two different approaches to estimate the admissible set of controls $\mathcal{U}_{n}$ in equation (3):

• •

Probability estimation: Given a state $\mathbf{X}_{n}=\mathbf{x}$ and $u\in\mathcal{W}$, we estimate, via simulation, the probability of violating the constraint

[TABLE]

It follows that $u\in\hat{\mathcal{U}}_{n}(\mathbf{x})\;\Leftrightarrow\;\hat{p}_{n}(\mathbf{x},u)<p$. Particularly, to compute $\hat{p}_{n}(\mathbf{x},u)$ we consider Gaussian process smoothing of empirical probabilities, logistic regression and parametric density fitting.

• •

Quantile estimation: We approximate the quantile $q_{n}(\mathbf{x},u)$ of $G_{n}(\mathbf{x},u)$ via empirical ranking, support vector machines and quantile regression methods. The admissible sets $\mathcal{U}_{n}(\mathbf{x})$ and $\mathcal{X}^{a}_{n}(u)$ are then defined as:

[TABLE]

To implement all of the above techniques we use Monte Carlo simulation, specifying first the simulation design and then sampling (independently across draws) the $G$’s or $Y$’s to be used as training data. We work in a flexible framework where samples of $G_{n}(\mathbf{x},u)$ are generated in batches of $M_{b}$ simulations from each design site $\{\mathbf{x}^{i},u^{i}\}_{i=1}^{M_{a}}$. The case of $M_{b}=1$ corresponds to a classical regression approach, while large $M_{b}\gg 1$ can be interpreted as nested Monte Carlo averaging along $M_{b}$ inner samples.

Remark 4.1**.**

In section 3.2, we parameterized the elements of the approximation space $\mathcal{H}_{n}^{c}$ for estimation of the continuation value function $\hat{\mathcal{C}}(\cdot,\cdot)$ via vectors $\bm{\alpha}$ i.e. $f_{n}^{c}(\mathbf{x},u)\equiv f_{n}^{c}(\mathbf{x},u;\bm{\alpha})$, cf. (24) and (25). To distinguish, in the following sections we use $\bm{\beta}$ to parameterize the approximators in $\mathcal{H}_{n}^{a}$ for estimating the admissible set: $f_{n}^{a}(\mathbf{x},u)\equiv f_{n}^{a}(\mathbf{x},u;\bm{\beta})$ in eq. (17). The meaning and dimension of $\bm{\beta}$ will vary from method to method.

4.1 Probability estimation

4.1.1 Interpolated nested Monte Carlo (INMC)

Recall the NMC method from Section 2.1 where we select $M_{a}$ design sites of state-action pairs and simulate multiple paths from each site to locally assess the probability of $G_{n}(\mathbf{x},u)>0$ (in what follows, we suppress in the notation the dependence on $n$). Specifically, for each design site $(\mathbf{x}^{i},u^{i}),\ i=1,\ldots,M_{a}$, we simulate $M_{b}$ batched samples from the distribution $G(\mathbf{x}^{i},u^{i})$ as $\{g^{b}(\mathbf{x}^{i},u^{i})\}_{b=1}^{M_{b}}$. The unbiased point estimator of $p(\mathbf{x}^{i},u^{i})$ is:

[TABLE]

Since (33) only yields $M_{b}$ local estimates $\bar{p}(\mathbf{x}^{i},u^{i})$, for Algorithm 1 we have to extend them to an arbitrary state-action $(\mathbf{x},u)\mapsto\hat{p}_{\text{INMC}}(\mathbf{x},u)$. This is achieved by interpolating $\bar{p}(\mathbf{x}^{i},u^{i})$’s, e.g. linearly. The admissible set with confidence level $\rho$ becomes:

[TABLE]

However, especially for $M_{b}$ small, interpolation performs poorly because the underlying point estimates $\bar{p}(\mathbf{x}^{i},u^{i})$ are noisy. Therefore, smoothing should be applied via a statistical regression model. Regression borrows information cross-sectionally to mitigate the estimation noise, reducing the variance of $\bar{p}$.

4.1.2 Gaussian process regression (GPR)

GPR is a flexible non-parametric regression method that views the map $(\mathbf{x},u)\rightarrow p(\mathbf{x},u)$ as a realization of a Gaussian random field so that any finite collection of $\{p(\mathbf{x},u)\},(\mathbf{x},u)\in\mathcal{X}\times\mathcal{W}$ is multivariate Gaussian. For any $n$ design sites $\{(\mathbf{x}^{i},u^{i})\}_{i=1}^{n}$, GPR posits that

[TABLE]

with mean vector $\vv{m}_{n}:=[m(\mathbf{x}^{1},u^{1};\bm{\beta}),\ldots,m(\mathbf{x}^{n},u^{n};\bm{\beta})]$ and $n\times n$ covariance matrix $\mathbf{K}_{n}$ comprised of $\kappa(\mathbf{x}^{i},u^{i},\mathbf{x}^{i^{\prime}},u^{i^{\prime}};\bm{\beta}),\text{ for }1\leq i,i^{\prime}\leq n$. The vector $\bm{\beta}$ represents all the hyperparameters for this model.

Given the training dataset $\{(x^{i},u^{i}),\bar{p}^{i}\}_{i=1}^{M_{a}}$ (where $\bar{p}^{i}$ is a shorthand for $\bar{p}(\mathbf{x}^{i},u^{i})$), GPR infers the posterior of $p(\cdot,\cdot)$ by assuming an observation model of the form $\bar{p}(\mathbf{x},u)=p(\mathbf{x},u)+\epsilon$ with a Gaussian noise term $\epsilon\sim\mathcal{N}(0,\sigma^{2}_{\epsilon})$. Conditioning equations for multivariate normal vectors imply that the posterior predictive distribution $p(\mathbf{x},u)|\{(x^{i},u^{i}),\bar{p}^{i}\}_{i=1}^{M_{a}}$ at any arbitrary site $(\mathbf{x},u)$ is also Gaussian with the posterior mean $\hat{p}_{\text{GPR}}(\mathbf{x},u)$ that is the proposed estimator of $p(\mathbf{x},u)$:

[TABLE]

and $\mathbf{K}$ is $M_{a}\times M_{a}$ covariance matrix described through the kernel function $\kappa(\cdot,\cdot;\beta)$.

The mean function is often assumed to be constant $m(\mathbf{x},u;\bm{\beta})=\beta_{0}$ or described using a linear model $m(\mathbf{x},u;\bm{\beta})=\sum_{k=1}^{K}\beta_{k}\phi(\mathbf{x}^{i},u^{i})$ with $\phi(\cdot,\cdot)$ representing a polynomial basis. A popular choice for the kernel $\kappa(\cdot,\cdot,\cdot,\cdot)$ is squared exponential (see equation (36)) with $\{\{\beta_{\mathrm{len},k}\}_{k=1}^{d},\beta_{\mathrm{len},u}\}$ termed the lengthscales and $\sigma_{p}^{2}$ the process variance of $p(\cdot,\cdot)$:

[TABLE]

The hyperparameters $\bm{\beta}:=(\{\beta_{k}\}_{k=1}^{K},\{\beta_{\mathrm{len},k}\}_{k=1}^{d},\beta_{\mathrm{len},u},\sigma^{2}_{p},\sigma^{2}_{\epsilon})$ represent different attributes of $\hat{p}_{\text{GPR}}(\cdot,\cdot)$. The parameters $\{\beta_{k}\}_{k=1}^{K}$ determine the trend of $\hat{p}(\cdot,\cdot)$ over the input domain. The lengthscales $\{\beta_{\mathrm{len},k}\}_{k=1}^{d}$ and $\beta_{\mathrm{len},u}$ determine spatial smoothness or how quickly the function changes. Small lengthscales make $\hat{p}$ to be “bumpy”, while large lengthscales make $\hat{p}$ close to linear. The process variance $\sigma^{2}_{p}$ determines the amplitude of fluctuations in $\hat{p}_{\text{GPR}}$ and $\sigma^{2}_{\epsilon}$ represents the sampling noise variance. These hyperparameters are estimated by maximizing the log-likelihood function based on the dataset $\{(x^{i},u^{i}),\bar{p}^{i}\}_{i=1}^{M_{a}}$. Besides squared exponential kernel (Equation 36), other popular kernels include Matérn-3/2 and Matérn-5/2 [35].**

A conservative estimate $\hat{p}_{\text{GPR}}^{(\rho)}(\mathbf{x},u)$ at confidence level $\rho$ is obtained by explicitly incorporating the (estimated) standard error of $\bar{p}(\mathbf{x}^{i},u^{i})$ into the GPR smoothing. Namely, we adjust the training dataset to $\{(x^{i},u^{i}),\bar{p}^{i}_{\rho}\}_{i=1}^{M_{a}}$, where $\bar{p}^{i}_{\rho}:=\bar{p}(\mathbf{x}^{i},u^{i})+z_{\rho}\sqrt{\frac{\bar{p}(\mathbf{x}^{i},u^{i})(1-\bar{p}(\mathbf{x}^{i},u^{i}))}{M_{b}}}$. The resulting $\hat{p}_{\text{GPR}}^{(\rho)}(\mathbf{x},u)$ is the counterpart of (34) using $\{(x^{i},u^{i}),\bar{p}^{i}_{\rho}\}_{i=1}^{M_{a}}$.

In Figure 2(b) we present the dataset $\{L^{i},I^{i},0,\bar{p}^{i}\}_{i=1}^{M_{a}}$ (background colormap) for the microgrid case study. The thick red line indicates the contour $\{\hat{p}_{GPR}=5\%\}$, dividing the state space $\mathcal{X}$ for $u=0$ into admissible $\mathcal{X}^{a}(0)$ (left of red line) and inadmissible region $(\mathcal{X}^{a}(0))^{c}$ (right of red line).

4.1.3 Logistic regression (LR)

In the previous section, we created local batches to estimate ${p}(\mathbf{x}^{i},u^{i})$ pointwise and then regressed these estimates to build a global approximator. A classical alternative is to learn the probability of $G(x,u)>0$ using a logistic regression model. This setup uses a single sample $g(\mathbf{x}^{i},u^{i})$ from $G(\mathbf{x}^{i},u^{i})$ from each design site $(\mathbf{x}^{i},u^{i})$ and transforms it to a binary response $y^{i}=\mathbbm{1}_{g(\mathbf{x}^{i},u^{i})>0}$. The probability $\hat{p}(\mathbf{x},u)$ is then modeled as a generalized linear model with a logit link function

[TABLE]

The basis functions $\phi(\mathbf{x},u)$ could be polynomials, e.g. quadratic or cubic in coordinates of $(\mathbf{x},u)$. The regression coefficients ${\bm{\beta}}$ are fitted using the dataset $\{\mathbf{x}^{i},u^{i},\mathbf{y}^{i}\}_{i=1}^{M_{s}}$, as the solution to

[TABLE]

We may again create a more conservative estimate $\hat{\mathcal{U}}^{(\rho)}_{LR}(\mathbf{x})$ of $\hat{\mathcal{U}}_{LR}(\mathbf{x})$ at confidence level $\rho$ by utilizing the standard error for $\hat{p}_{LR}$ using the Delta method **[40]**:

[TABLE]

In Figure 2(a), we present the original realizations $y^{i}\in\{0,1\}$ (in blue) for a design in the input subspace $(L,I,u=0)$ of the microgrid case study. The figure indicates the resulting logistic regression fit $\hat{p}_{LR}(L,I,0)$ at levels $1\%$, $5\%$ and $10\%$ (i.e. contour lines of $\hat{p}_{LR}(\hat{\bm{\beta}})\in\{0.01,0.05,0.1\}$). The admissibility set for $u=0$, $\mathcal{X}_{n}^{a}(0)$ is the region to the left of the thick red contour.

Remark 4.2**.**

*Similar to Section 4.1, we can simulate batched samples from each design site for the logistic regression, leading to “binomial” observation likelihood instead of (38). *

Remark 4.3**.**

*A non-parametric variant of equation (37) is kernel logistic regression, where the basis functions are $\phi_{j}(\mathbf{x},u)=\kappa(\mathbf{x},u,\mathbf{x}^{j},u^{j})$ for a kernel function $\kappa$ centered at $(\mathbf{x}^{j},u^{j})$. One common choice is radial basis functions (RBF) where $\kappa(\mathbf{x},u,\mathbf{x}^{j},u^{j})=\exp{(-\gamma_{1}\|\mathbf{x}-\mathbf{x}^{j}\|^{2}_{2}-\gamma_{2}\|u-u^{j}\|^{2}_{2})}.$ RBF can be interpreted as the squared-exponential kernel for a logistic Gaussian Process model, with a fixed bandwidth parameter $\gamma_{i}$. In contrast, in GPR the bandwidths are estimated through MLE. *

4.1.4 Parametric density fitting (PF)

This approach fits the distribution $G(\mathbf{x},u)$, and then analytically infers the probability \mathbb{P}\Big{(}G(\mathbf{x},u)>0\Big{)} from the corresponding cumulative distribution function. This is done by proposing a parametric family $\{f(\cdot;\Theta)\}$ of densities, fitting the underlying parameters $\Theta$ based on an empirical sample from $G$ and then evaluating the resulting analytical probability $\bar{p}_{PF}(\mathbf{x},u):=\int_{0}^{\infty}f_{G(\mathbf{x},u)}(z|\hat{\Theta}(\mathbf{x},u))dz$. PF yields a “universal” solution across a range of constraint levels $p$.

At a design site $(\mathbf{x},u)$, the probability $p(\mathbf{x},u)$ is estimated in a two-step procedure: first estimated locally over a design $\mathcal{D}^{a}=\{\mathbf{x}^{i},u^{i}\}$ and then regressed over the full input domain $\mathcal{X}\times\mathcal{W}$. For the first step, we apply nested Monte Carlo to generate a collection of realized $\{g^{b}(\mathbf{x}^{i},u^{i})\}_{b=1}^{M_{b}}$ that is used to construct a parametric density via the maximum likelihood estimate:

[TABLE]

In the second step, we evaluate $\tilde{p}_{PF}(\mathbf{x}^{i},u^{i}):=\int_{0}^{\infty}f_{G}(z|\hat{\Theta}(\mathbf{x}^{i},u^{i}))$ and extend it to the full domain $\mathcal{X}\times\mathcal{W}$ based on the computed $\{\mathbf{x}^{i},u^{i},\tilde{p}_{PF}(\mathbf{x}^{i},u^{i})\}_{i=1}^{M_{a}}$ using $\mathcal{L}_{2}$ projection:

[TABLE]

where $\mathcal{M}_{T}$ is an approximation space chosen for regression. The admissible set $\mathcal{U}(\mathbf{x})$ is estimated as:

[TABLE]

A transformation of the distribution $G(\mathbf{x},u)$ might be important for above distribution fitting. For example, in the context of microgrid, in Section 2.2, G=\mathcal{L}\big{(}\sup_{s\in[t_{n},t_{n+1})}S(s)\big{)} has a point mass at [math] and thus, any continuous distribution will lead to poor statistical estimation. Using a transformation that preserves the probability of the target event,

[TABLE]

we work with G^{\prime}(L_{n},I_{n},u_{n}):=\mathcal{L}\big{(}\sup_{s\in[t_{n},t_{n+1})}[L(s)-u_{n}-\frac{I(s)}{\delta s}\wedge B_{\max}]\big{)}. In Figure 2(c) we present the empirical and estimated probability $z\mapsto\mathbb{P}(G^{\prime}(L_{n},I_{n},u_{n})>z)$ when $L_{n}=5.5,I_{n}=1.48$ and $u_{n}\in\{0,1\}$ for the microgrid example. We model the distribution $G^{\prime}$ using a truncated normal distribution, $\mathbb{P}(G^{\prime}\leq g)=\Phi(\frac{g-\theta_{2}}{\theta_{3}})\mathbbm{1}_{g\geq\theta_{1}}$, with parameters $\Theta=(\theta_{1},\theta_{2},\theta_{3})$ representing the location of censoring, the mean and the standard deviation respectively. At $L_{n}=5.5,I_{n}=1.48,u_{n}=1.0$ and inner simulation budget $M_{b}=100$, the estimated parameters $(\hat{\theta}_{1},\hat{\theta}_{2},\hat{\theta}_{3})=(-1.5,-1.12,0.53)$ result in probability $\tilde{p}_{PF}(5.5,1.48,1.0)=0.016$. The corresponding probability after $\mathcal{L}_{2}$ projection (equation (40)) is $\hat{p}_{PF}(5.5,1.48,1.0)=0.017$. Thus at $p=0.05$, the control $u=1.0\in\hat{\mathcal{U}}_{n}$ is admissible. However, at $u_{n}=0$, $(\hat{\theta}_{1},\hat{\theta}_{2},\hat{\theta}_{3})=(-0.5,-0.12,0.55)$, $\tilde{p}_{PF}(5.5,1.48,0.0)=0.414$ and $\hat{p}_{PF}(5.5,1.48,0.0)=0.429$, thus the control $u=0\notin\hat{\mathcal{U}}_{n}$ is inadmissible.

4.2 Quantile estimation

In this section we consider methods for modeling and estimating $q(\mathbf{x}^{i},u^{i})$, the $(1-p)$-th quantile of the distribution $G(\mathbf{x}^{i},u^{i})$. Admissibility corresponds to the quantile being negative.

4.2.1 Empirical percentiles (EP)

As before, we start by choosing $M_{a}$ design sites of state-action pairs and generate batched samples $\{g^{b}(\mathbf{x}^{i},u^{i})\}_{b=1}^{M_{b}}$ from each design site $(\mathbf{x}^{i},u^{i})$. The empirical estimate of $q(\mathbf{x}^{i},u^{i})$ is simply the $(1-p)^{th}$ percentile of the realized $\{g^{b}\}_{b=1}^{M_{b}}$ (which requires $M_{b}>p^{-1}$):

[TABLE]

Similar to previous methods, we extend to arbitrary $(\mathbf{x},u)\mapsto\hat{q}(\mathbf{x},u)$ using regression on the dataset $\{\mathbf{x}^{i},u^{i},\bar{q}(\mathbf{x}^{i},u^{i})\}_{i=1}^{M_{a}}$ and an approximation space $\mathcal{M}_{q}$. The set of admissible controls for $\mathbf{x}$ is: \hat{\mathcal{U}}_{EP}(\mathbf{x}):=\Big{\{}u:\hat{q}(\mathbf{x},u)\leq 0\Big{\}}. In Figure 2(e) we show the estimated $\hat{q}(\cdot,\cdot,0)$ indicated via the background colormap. The thick red line indicates the zero-contour $\hat{q}=0$, so that the admissibility set for $u=0$, $\mathcal{X}_{n}^{a}(0)$, is the region to the left of the contour.

Remark 4.4**.**

*This approach is similar to the INMC approach discussed in Section 4.1, however, here we model the quantile rather than the probability of exceeding zero. Furthermore, we can use the regression standard error of $\hat{q}(\cdot,\cdot)$ to construct a more conservative estimate of the admissible set $\mathcal{U}_{EP}(\mathbf{x})$. *

A popular alternative to adjusting $\bar{q}$’s via regression standard errors is to replace the empirical percentile with the empirical conditional tail expectation (CTE):

[TABLE]

and then running a regression on the training set $(\mathbf{x}^{i},u^{i},\overline{\text{CTE}}(\mathbf{x}^{i},u^{i}))$ to obtain the fitted CTE surface $\widehat{\text{CTE}}(\mathbf{x},u)$ and finally \hat{\mathcal{U}}_{CTE}(\mathbf{x}):=\big{\{}u:\widehat{\text{CTE}}(\mathbf{x},u)\leq 0\big{\}}. This idea is similar to regularizing the Value-at-Risk estimation with the Conditional VaR.

4.2.2 Support Vector Machines (SVM)

For a fixed control $u$, finding the admissible set $\mathcal{X}^{a}_{n}(u)$ in (9) can be interpreted as classifying each input $\mathbf{x}$ as being in $\mathcal{X}^{a}_{n}(u)$ or not. Therefore, we consider the use of classification techniques, specifically support vector machines (SVM). This approach does not estimate the $(1-p)$-quantile $q(\mathbf{x},u)$, but rather its [math]-level set with respect to $(\mathbf{x},u)$. The starting point is to use the nested Monte Carlo simulations to compute $\bar{p}(\mathbf{x}^{i},u^{i})$ with much smaller batch size $M_{b}$ compared to Section 4.1. Next, we construct a binary classification objective with a training dataset $\{\mathbf{x}^{i},u^{i},y^{i}\}_{i=1}^{M_{a}}$ where the $\pm 1$-labels are

[TABLE]

The boundary separating the two classes is evaluated by solving the optimization problem:

[TABLE]

where \phi(\mathbf{x},u)=\big{[}\phi_{1}(\mathbf{x},u),\phi_{2}(\mathbf{x},u),\ldots,\phi_{K}(\mathbf{x},u)\big{]}^{T} are the $K$ basis functions and $C$ is the penalty parameter. We estimate the set of admissible controls corresponding to $\mathbf{x}$ as:

[TABLE]

Figure 2(d) displays the estimated $\hat{\mathcal{X}}^{a}_{n}(u)$ and the corresponding dataset $(L^{i},I^{i},0,y^{i})$ ($u=0$* is fixed). The region where $u=0$ is admissible is to the left of the (thick red) decision boundary.*

Remark 4.5**.**

A conservative estimate $\hat{\mathcal{U}}^{(\rho)}_{SVM}$ is obtained by biasing the decision boundary to the left by re-labeling the training points in (42) via:

[TABLE]

4.2.3 Quantile Regression (QR)

QR directly constructs a parametric model for $q(\mathbf{x},u)$:

[TABLE]

To estimate the coefficients $\bm{\beta}\in\mathbb{R}^{K}$, we use the dataset $\{\mathbf{x}^{i},u^{i},g^{i}\}_{i=1}^{M_{a}}$ (where $g^{i}$ is a sample from the distribution $G(\mathbf{x}^{i},u^{i})$) to maximize the negative log likelihood:

[TABLE]

where $\mathcal{L}^{(p)}(y)=y(p-1_{\{y<0\}})=py_{+}+(1-p)y_{-}$. As for the parametric density fitting, a transformation of $G(\mathbf{x},u)$ might be beneficial when applying quantile regression. Figure 2(f) presents the dataset $\{L^{i},I^{i},0,g^{i}\}_{i=1}^{M_{a}}$ (background colormap) and the estimated admissible set $\hat{\mathcal{X}}^{a}(0)$ which is the region to the left of the contour $\{\hat{q}_{QR}(L,I)=0\}$ (thick red line).

Relying on the Delta method again to compute the variance of the estimated quantile $\hat{q}(\mathbf{x},u;\hat{\bm{\beta}})$ as $\phi(\mathbf{x},u)^{\prime}Var(\hat{\bm{\beta}})\phi(\mathbf{x},u)$, the admissible set at $\mathbf{x}$ at confidence level $\rho$ is:

[TABLE]

5 Case Studies

Recall the problem introduced in Section 2.2 where we control the operation of a diesel generator in order to match demand at minimal cost while maintaining the probability of blackout between each decision epoch below a given threshold $p$. In this section, we discuss two variants of such microgrid control. In the first example, we assume a time-homogeneous net-demand process which reduces the problem of estimating admissible set to a pre-processing step. In the second example, we use time-dependent net demand process calibrated to data obtained from a microgrid in Huatacondo, Chile. Time-inhomogeneity requires to estimate the admissible set at every step. The microgrid features a perfectly efficient battery, so that the respective power output at $t_{n_{k}}$ (recall $t_{n_{k}}$ is a generic time instance on the finely discretized time grid) is given by:

[TABLE]

Table 1 lists other microgrid parameters, i.e. capacity of the battery $I_{\max}$, maximum charging rate $B_{\min}$, maximum discharging rate $B_{\max}$ and diesel switching cost $\mathcal{K}$.

5.1 Implementation details

Numerical Gold Standard:* In the absence of analytic benchmark, we use a high-budget gold standard to compare the output from the models discussed in Section 4. For each fixed time-step $t_{n}$ we discretize the domain $\mathcal{X}=(L,I)$ into $10,000$ design sites over a grid of $100\times 100$. For each design site $(L^{i},I^{j}),\ i,j\in\{1,\ldots,100\}$ and $u^{k}\in 0\cup\{1=u_{1},\ldots,u_{101}=10\}$, we evaluate $\hat{p}(L^{i},I^{j},u^{k})$ using (33) with batch size $M_{b}=10,000$. Thus, the total simulation budget is $100\times 100\times 102\times 10000\approx 10^{10}$. We then evaluate the local minimal admissible control*

[TABLE]

To evaluate $u^{\min}_{n}(L,I)$ at new sites we employ linear interpolation on the dataset $\{L^{i},I^{j},u^{\min}_{n}(L^{i},I^{j})\}_{i,j=1}^{100}$.

To estimate the continuation function, we use the piecewise continuous approximation of Section 3.2 with $M_{I}=15,M_{u}=15$ and $M_{L}=1500$ sites in $L$. The design $\mathcal{D}^{c}$ is constructed as the Cartesian product $\{L^{1},L^{2},\ldots,L^{M_{L}}\}\times\{I^{0},I^{1},\ldots,I^{M_{I}}\}\times\{u^{0},u^{1},\ldots,u^{M_{u}}\}$, where $L^{i},i=1,\ldots,M_{L}$ are sampled uniformly from the interval $[L_{\min},L_{\max}]=[-8,8]$. The inventory $\{I^{0},I^{1},\ldots,I^{M_{I}}\}$ and control $\{u^{1},\ldots,u^{M_{u}}\}$ discretizations are equispaced in $[0,I_{\max}]$ and $[\underline{u},\overline{u}]$ respectively.

Low budget policies:* We approximate the continuation value function $\mathcal{C}$ using a piecewise continuous approximation with degree-3 in $L$ combined with interpolation in other dimensions. The design $\mathcal{D}^{c}$ is the same as for the numerical gold standard with discretization levels $M_{L}=1000,M_{I}=10,M_{u}=10$.*

We approximate the admissible set $\mathcal{U}$ using the methods described in Section 4 and compare the performance of each method by using a fixed set of $M^{\prime}=20,000$ out-of-sample simulations. To address the discontinuity in $\mathcal{W}=0\cup[\underline{u},\overline{u}]$, we implement two separate statistical models to learn $\mathcal{U}_{n}(\cdot)$. As an example, with logistic regression of Section 4.1.3 we estimate two sets of parameters in equation (37): the first one uses one-step paths generated with $u=0$ and a two-dimensional regression of $y^{i,(1)}$ against $(L^{i},I^{i})$. The second one uses design sites in the three-dimensional space $(L,I,u)$ where $u\in[1,10]$ and a 3-D regression of $y^{i,(2)}$ against $(L^{i},I^{i},u^{i})$. For both regressions we choose a Sobol sequence space-filling experimental design: $\mathcal{D}^{a}=(L^{i},I^{i})\in[-2,8]\times[0,10]$ when $u=0$, and $\mathcal{D}^{a}=(L^{i},I^{i},u^{i})\in[-2,8]\times[0,10]\times[1,10]$ otherwise. The control space $[1,10]$ is discretized into 51 levels.

Additional parameters used for each method are specified in Table 2. We found that Matern-3/2 kernels work better than (36) for smoothing $\bar{p}(L,I,u)$ (GPR) and $\tilde{p}(L,I,u)$ (PF) because the respective input-output maps feature steep transitions as a function of $(L,i,u)$. It is known that “rougher” kernels are better suited for such learning tasks compared to the $C^{\infty}$-smooth squared exponential kernel (36) by allowing the fitted $\hat{p}$ to have more “wiggle room”. On the other hand, in the context of EP and CTE the input observations of $\hat{q}(L,I,u)$ and $\overline{\text{CTE}}(L,I,u)$ are quite smooth in $(L,i,u)$ and both GP kernel families perform equally well. The algorithms are implemented in python 2.7. We used “GaussianProcessRegressor” and “SVM.SVC” functions from sklearn library for GPR and SVM respectively. For LR and QR we used “Logit” and “quantile_regression” functions from statsmodels library.

5.2 Example 1: Microgrid with Stationary Net Demand

In this subsection, we assume time-homogeneous Ornstein-Uhlenbeck dynamics of the net demand process

[TABLE]

where $(B(t))$ is a standard Brownian motion. This scenario reduces the complexity of learning the probability constraints since we need to estimate the admissible set $\mathcal{U}_{0}(\cdot)$ only once as a pre-processing step before starting the approximate dynamic programming scheme for the continuation values. The simplified setting offers a good testbed to evaluate the performance of different admissible set estimation methods of Section 4; we show that the relative performance remains similar as we extend to more realistic dynamics in Section 5.3. For this example, we assume the mean reversion parameter $\lambda=0.5$ and volatility $\sigma=2$.

Figure 3(a) plots the resulting costs $\hat{V}_{0}(0,5)$ versus the frequency of inadmissible decisions $w_{freq}$ for different methods of Section 4. We show the results both for $p=0.05$ (dark blue), and $p=0.01$ (light grey) and benchmark both cases against the numerical gold standard. Since the probabilistic constraints form the crux of the problem, we require schemes to maintain $\hat{u}\in\mathcal{U}_{n}$ as much as possible, i.e., $w_{freq}\approx 0$. At $p=5\%$, we observe 0.09% , 0.54% and 1.36% frequency of inadmissible decisions with logistic regression (LR), Gaussian process regression (GPR) and parametric density fitting (PF), respectively. Such accuracy might be deemed acceptable. However $w_{freq}$ is much worse (as high as 8.4% with EP) for the other methods. While all the methods are a priori consistent, admissible set estimation via probability-based methods clearly seems to outperform quantile-based ones. Our experiments suggest that at low simulation budget, estimators of ${p}(\mathbf{x},u)$ have significantly lower bias compared to estimators of ${q}(\mathbf{x},u)$, thus partially explaining the difference. For a more stringent threshold $p=1\%$, we find the cost of all the methods to increase, without significant difference in the frequency of inadmissible decisions $w_{freq}$. Indeed, Figure 3 illustrates the trade-off between lower costs and lower $w_{freq}$ (i.e. more conservative estimate of the constraints).

Table 3 expands Figure 3 by also reporting the corresponding $\tilde{\mathcal{T}}$ statistic, the average inadmissibility margin $w_{avm}$ and realized frequency of violations (i.e. blackouts) $w_{rlzd}$ defined as:

[TABLE]

We find the realized frequency of violations $w_{rlzd}$ to be lowest for LR, GPR and PF. The average inadmissibility margin $w_{avm}$ is also lowest for GPR and PF (the large value of $w_{avm}$ for LR is attained in very small region as evident from $w_{freq}\approx 0$). The $\tilde{\mathcal{T}}$ statistic is negative for LR, GPR and PF and positive for the rest, meaning that all other methods fail to statistically respect the probability constraints when binding. Due to small frequency of inadmissible decisions $w_{freq}$, cost $\hat{V}_{0}(0,5)$ similar to the numerical gold standard and negative test statistic $\tilde{\mathcal{T}}$, we recommend LR, GP and PF methods for the problem at hand.

Next, we test the sensitivity of the cost in terms of the probability threshold $p$ (employing logistic regression $\hat{\mathcal{U}}_{LR}$) in Figure 3(b). Increasing $p$ decreases $V$ as the set of admissible controls $\mathcal{U}$ monotonically increases in $p$. For example, any admissible control at $p=1\%$ threshold is also feasible for $p>1\%$, thus the respective cost at $1\%$ will be at least as much as at, say, $10\%$ threshold.

As previously discussed, the constraint is binding for only approximately $10\%$ of time-steps. In fact, that probability varies across the methods since the estimate of $\hat{\mathcal{U}}$ affects the choice of $\hat{u}_{n}$ and ultimately the distribution of $\hat{X}_{n}$. Intuitively, the realized system states are driven by the estimates of the probabilistic constraints. Typically, more conservative estimates of $\mathcal{U}$ will push $\hat{X}_{0:N}$ away from the “risky” regions. This is also confirmed in Figure 3 where as $p\rightarrow 1$, $w_{rlzd}\rightarrow 20\%=w_{bind}$ while in Table 3 $w_{bind}\simeq 10\%$.

The metrics $w_{freq}$ (28), $w_{bind}$ (32), $w_{rlzd}$ (47) are closely linked. As the inadmissible decisions can occur only when the constraint is binding, $u^{\min}>0$, we expect $w_{freq}\leq w_{bind}$ and $w_{freq}\approx w_{bind}$ for a method with a bias in overestimating the admissible set (e.g. $\mathcal{X}^{a,EP}(u)\supset\mathcal{X}^{a,GS}(u)\ \forall u\in\mathcal{W}$). The realized violations (blackouts) $w_{rlzd}$ can be represented as a sum of three:

[TABLE]

where the weights $p_{1},p_{2},p_{3}$ depend on the distribution of the controlled trajectories. The first term represents the instances when the constraint is binding but the controller chooses an inadmissible control (i.e. mis-estimates $\hat{\mathcal{U}}$). The second term represents instances when the constraint is binding and correctly estimated, but due to random shocks violations take place (with a conditional frequency below the specified $p=0.05$). The last term represents instances when the constraint is not binding but some violations still occur with the intrinsic conditional frequency strictly less than $p$. Note that due to $w_{bind}\ll 1$, most of the violations are of the latter type, i.e. take place when $u^{*}=0$ and the conditional violation probability is below $p$. We illustrate these scenarios in Figure 3(c) using the LR model. Thus, the first term counts the instances when violations occur at the same time as controller makes an inadmissible decision (circle encircling triangle), the second term counts the triangles when $I\approx 0$, and the third term the triangles in the grey region where the constraint is not binding (violations when $u^{\min}=0$).

Although we observed poor performance of quantile based methods, asymptotically (with respect to the simulation budget) we expect them to perform similar to the probability based methods. As an example, in Appendix A Table 6, we present the performance of SVM for thresholds $p=5\%$ and $p=1\%$ with increasing budget. For $p=5\%$ and by increasing the simulation budget from $10^{5}$ to $10^{8}$, we find the frequency of inadmissible decisions $w_{freq}$ to drop from $5.93\%$ to $1.5\%$, average inadmissibility amount $w_{avm}$ from $0.78$ kW to $0.27$ kW, frequency of realized blackouts $w_{rlzd}$ from $2.80\%$ to $0.30\%$ and the test statistic which rejected the method at $10^{5}$ simulation budget ($\mathcal{T}\gg 0$) suggests to accept it ($\mathcal{T}\ll 0$) at $10^{8}$ simulation budget. We observe similar behavior at $p=1\%.$ The main challenge with quantile-based methods is the underlying bias in learning $q(\mathbf{x},u)$. This bias is known to converge to zero slowly, necessitating a relatively large $M_{b}$. A glimpse of this can be observed in Appendix A where SVM does not perform adequately all the way up to total budget of $10^{8}$. Increasing $M_{b}$ with $M_{a}$ fixed offers only a limited improvement that tapers off quickly, as was observed previously in[26]. Keeping a constant budget, a high $M_{b}$ forces a low $M_{a}$ which offsets these performance gains as the regression is not able to properly explore the space.

There are many techniques to improve $\bar{q}$, for example variance reduction tools. We have experimented with variance reduction for the SVM method by using antithetic variables. Although the results were marginally better, they were not meaningfully different from those reported in Appendix A. A completely different way to improve EP and SVM could be to judiciously choose the simulation design $\mathcal{D}^{a}$. Paper [26] explores various approaches to build simulation designs for estimating the continuation value function and documents their strong impact on performance. We anticipate that a similar strategy for $\mathcal{D}^{a}$ may improve convergence.**

Conservative estimators for $\mathcal{U}$.* Algorithms for SCPC are expected to respect the probabilistic constraints, so that it is critical to minimize the occurrence of inadmissible decisions. As discussed in Section 4, one way to raise the statistical guarantee for admissibility of $\hat{\mathcal{U}}$ is by adding a margin of error $\xi(\mathbf{x},u)$. This yields a more conservative (i.e. smaller) $\hat{\mathcal{U}}$ and hence lowers $w_{freq}$. In Table 4 we examine three scenarios for $\xi(\mathbf{x},u)$ sorted from least to most conservative (in all cases we maintain $p=5\%$):*

• •

Scenario 1: unadjusted $\xi=0\%$ (same as Table 3);

• •

Scenario 2: $\xi^{(\rho)}(\mathbf{x},u)$ at 95% confidence level, $z_{\rho}=1.96$;

• •

Scenario 3: fixed $\xi=4\%$, which is equivalent to lowering the violation threshold to $p-\xi=1\%$.

Table 4 confirms the intuition that the frequency of inadmissible decisions $w_{freq}$ should be decreasing from scenario 1 to 3. This is illustrated in Figure 4 that shows how the minimum admissible control is affected by $\xi(\mathbf{x},u)$. Although adding a margin of error does lower $w_{freq}$, this mechanism does not really alter the relative performance of the different methods. Thus, for all three scenarios, we find SVM, CTE and EP to give unreliable estimates of $\mathcal{U}$ (since $\tilde{\mathcal{T}}\gg 0$). An exception is QR which yields high $w_{freq}$ for $\xi=0$ but does become acceptable ($\tilde{\mathcal{T}}<0$) in scenario 3. In contrast, LR, GPR and PF perform well throughout. Table 7 in Appendix B provides additional results as we vary the confidence level to $\rho=90\%,99\%\text{ and }99.95\%$, with the same general conclusions. (Observe that a $w_{freq}$ very close to zero likely implies that $\hat{\mathcal{U}}\subset\mathcal{U}$ is strictly smaller and the controller is so conservative as to rule out some admissible actions.)

We generally expect the ultimate cost $\hat{V}_{0}(0,5)$ to increase as $\hat{\mathcal{U}}$ becomes more conservative, see the estimated $\hat{V}$’s across each row of Table 4. The increase in costs arises due to two factors: when the diesel generator is started sooner (due to $u=0$ becoming inadmissible as $\xi$ is raised), and the higher level of $\hat{u}$ once the diesel is ON. This can be seen in Figure 4 where in Scenarios 2 and 3 the controller switches the generator at a lower net demand and once the diesel is running picks a higher power output ($\hat{u}^{\min}(\cdot,I;p=5\%,\xi)-\hat{u}^{\min}(\cdot,I;p=5\%,\xi=0)>0$). We stress that the link between $\hat{\mathcal{U}}$ and $\hat{V}$ is complicated by the fact that as $\hat{\mathcal{U}}$ changes, so does the distribution of the controlled paths. So for example in Table 4 the cost for QR falls in Scenario 2, although it remains within two Monte Carlo standard errors.

Take-aways.* Our experiments demonstrate the following: (i) To accurately estimate admissible sets of the form (3) we recommend to use LR, GPR or PF which all model the underlying probability of violations $p(\mathbf{x},u)$. Although asymptotically equivalent, the approach of quantile estimation leads to poor estimates $\hat{\mathcal{U}}_{n}$ for practical budgets. (ii) Frequency of inadmissible decisions can be partly controlled by using conservative estimates $\hat{\mathcal{U}}^{\xi}_{n}$ at the expense of higher costs. However, even a conservative $\hat{\mathcal{U}}^{\xi}$ fails to make quantile-based methods acceptable, except for QR. (iii) For a new application, our suggested approach is to first evaluate the test statistic $\tilde{\mathcal{T}}$ at $\xi=0\%$ using one of the recommended methods. Depending on how close is $\tilde{\mathcal{T}}$ to zero, one can then adjust $\hat{\mathcal{U}}$’ via $\xi$ to improve the statistical guarantees on the frequency of inadmissible decisions $w_{freq}$.*

5.3 Example 2: Microgrid with seasonal demand

Unlike the previous example, where we assumed time-homogeneous net demand, in practice there is seasonality: during the day renewable generation is high and net demand is often negative; during morning/evening demand exceeds supply making $L(t)>0$. To incorporate this seasonality we use time-dependent Ornstein Uhlenbeck process (see **[17]** for a similar microgrid control problem):

[TABLE]

Here, $\lambda$ represents the speed of mean reversion towards the seasonal mean $\mu(t)$, while $\sigma(t)$ represents the time-varying volatility. Using Itô’s lemma and integration by parts one can prove that

[TABLE]

Thus,

[TABLE]

We calibrate $\mu(t)$ and $\sigma(t)$ in (67) using iterative methodology described in **[17]** and the data from a solar-powered microgrid in Huatacondo, Chile111https://microgrid-symposiums.org/microgrid-examples-and-demonstrations/huatacondo-microgrid/. Specifically, we compute the mean and variance of the residual demand over 24 hours at 15-minute intervals using data from Spring 2014, i.e. compute $\{\mu_{1},\mu_{2},\ldots,\mu_{96}\}$ and $\{\sigma_{1},\sigma_{2},\ldots,\sigma_{96}\}$. The estimated $\mu(t)$ can be seen in the left panel of Figure 5 that plots the empirical average of $L(t)$. As expected, during the day, i.e., $t\in[12,20]$ (noon-8:00 pm), the expected net-demand is negative ($\mu(t)<0$) while it is positive ($\mu(t)>0$) in the morning and during the night. The volatility $\sigma(t)$ is higher during the day due to the intermittent and unpredictable nature of solar irradiance. The mean reversion parameter was estimated to be $\lambda=0.3416$.

To visualize the interplay of the net demand, inventory and optimal control, the left panel of Figure 5 presents the average trajectories of the three processes over 48 hours. During the morning hours when the demand $L(t)$ is high and the battery is empty, the controller uses the diesel generator. During the day when the renewable output is high and $L(t)$ is negative, the controller switches off the diesel and the battery charges itself. However, the non-trivial region is when the average net-demand changes sign, either from positive to negative around noon or negative to positive in the evening. During the former time interval, the optimal control process is in $\{0,1\}$ (recall that minimum diesel output is $1$). Similarly, during the evening when the net demand becomes positive (as the renewable output declines), the controller quickly ramps up the diesel to match $L(t)\gg 0$. The right panel of Figure 5 repeats the average control and inventory curves, also showing their 2-standard deviation bands (in terms of the out-of-sample trajectories of $\hat{L}^{\hat{u}}_{0:T}$). As expected, the time periods around ramp-up or ramp-down of the diesel generator is when $\hat{u}_{n}$ experiences the greatest path-dependency and dispersion and differs most from the demand curve.

Comparing Table 5, which lists the estimated cost $\hat{V}_{0}(\mu(0),5)$ along with related statistics, with Figure 3 indicates that incorporating seasonal net-demand process does not change the relative order of performance between the methods. The cost goes up as the diesel generator has to be used throughout the mornings and the evenings to match demand.

As in the previous example, the performance of LR, GPR and PF almost matches the gold standard despite significantly lower simulation budget. In this setting the constraint is binding approximately 45% of the time (except for GPR and PF where it is 30% and 25% of the time). Frequency of inadmissible decisions $w_{freq}$ is 0.03% for LR, 1.17% for GPR, and 0.02% for PF. In contrast $w_{freq}$ is 43% for QR, 22% for EP, 43% for SVM and 22% for CTE, implying that all these schemes are highly unreliable for learning $\hat{\mathcal{U}}$. The average inadmissibility margin $w_{avm}$ is also significantly lower for GPR ($0.14$* kW) and PF ($0.26$ kW) compared to the rest of the methods. Here again we observe larger inadmissibility margin and very low frequency of inadmissible decisions for logistic regression. Similar behavior is also evident for the test statistic $\tilde{\mathcal{T}}$ and realized frequency of violations $w_{rlzd}$.*

To illustrate the typical behavior over a trajectory, Figure 6 plots the average control $Ave({\hat{u}_{n}}):=\frac{1}{M^{\prime}}\sum_{m=1}^{M^{\prime}}\hat{u}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})$ corresponding to different methods and the average minimum admissible control $Ave({u^{\min}_{n}}):=\frac{1}{M^{\prime}}\sum_{m^{\prime}=1}^{M^{\prime}}u^{\min}_{n}(\mathbf{x}_{n}^{\hat{u},m^{\prime}})$ computed using the gold standard. Notice that the latter is dependent upon the controlled trajectories $\mathbf{x}_{n}^{\hat{u}}$, resulting in different $Ave({u^{\min}_{n}})$ across methods. We expect $Ave({\hat{u}_{n}})$ above $Ave({u^{\min}_{n}})$ if a given method does not violate the constraint most of the time. This is true for LR and GPR, but SVM quite obviously fails, as the dashed line in the leftmost panel of Figure 6 is significantly higher than the solid line at numerous time steps. Furthermore, the conservative nature of GPR is reflected in the large difference between the average minimum admissible control and the average optimal control. This is also evident through $w_{bind}\approx 30\%$ for GPR compared to approximately 45% for the rest of the methods.

6 Conclusion

We developed a statistical learning framework to solve stochastic optimal control problems with local probabilistic constraints. The key objective of our algorithm is to efficiently estimate the set of admissible controls $\mathcal{U}(\cdot)$ and the continuation value function $\mathcal{C}(\cdot,\cdot)$ covering a general formulation of the state process dynamics and rewards. Since SCPC problems require estimating the admissible set repeatedly during the backward induction, we use regression based functional representation of $\mathbf{x}\mapsto\mathcal{U}(\mathbf{x})$. This perspective also provides a natural way of uncertainty quantification for admissibility, in particular offering conservative estimates that bring statistical guarantees regarding $\hat{\mathcal{U}}$. At the same time, our dynamic emulation algorithm allows parallel computation of $\mathcal{U}$ and $\mathcal{C}$ for additional computational efficiency.

Thanks to the plug-and-play functionality of the dynamic emulation algorithm, it was straightforward to test a large variety of schemes for learning $\mathcal{U}$. Our numerical results suggest that estimating probabilistic constraints via logistic regression, Gaussian process smoothing and parametric density fitting is more accurate than estimating the corresponding quantile (empirical ranking, SVM or quantile regression). A future line of research would be to additionally parametrize (e.g. using another GP model) the optimal control map $\mathbf{x}\mapsto\hat{u}_{n}(\mathbf{x})$ **[14]** which would speed-up the algorithm in the context of continuous action spaces. Another direction would be to extend the one-dimensional control framework described in this paper to a multi-dimensional setting. In reference to the microgrid example, multi-dimensional control will allow us to control the diesel output and the demand response or to control output from multiple dispatchable generators.

Appendix A Effect of Simulation Budget $(M_{b}\times M_{a})$

Appendix B Effect of Adaptive Margin of Error Level $\rho$