Input Hard Constrained Optimal Covariance Steering
Kazuhide Okamoto Panagiotis Tsiotras

	Improving the Scalability of a Prosumer Cooperative Game with K-Means Clustering

This work was supported in part by the Engineering and 
Physical Sciences Research Council under Grants EP/N03466X/1 and 
EP/S000887/1, and in part by the Oxford Martin Programme on Integrating 
Renewable Energy.
©2019 IEEE.  Personal use of this material is permitted.  Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. 
Digital Object Identifier 10.1109/PTC.2019.8810558.

    Liyang Han, Thomas Morstyn, Constance Crozier, Malcolm McCulloch
Department of Engineering Science 

University of Oxford

Oxford, United Kingdom 

{Liyang.Han, Thomas.Morstyn, Constance.Crozier, Malcolm.McCulloch}@eng.ox.ac.uk

    December 30, 2023
=============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

empty
	empty
	
	We address the optimal covariance steering (OCS) problem for stochastic discrete linear systems with additive  Gaussian noise under state chance constraints and input hard constraints. 
	Because the system state can be unbounded due to the unbounded noise, the state constraints are formulated as probabilistic (chance) constraints, i.e., the maximum probability of constraint violation is constrained.
	In contrast, because it is hard to interpret the appropriate control action when the control command violates the constraints, probabilistically formulating the control constraints are difficult, and deterministic hard constraints are preferable.
	In this work we introduce an OCS approach subject to simultaneous state chance constraints and input hard constraints and validate the approach using numerical simulations.
	



§ INTRODUCTION


The optimal covariance steering (OCS) problem is a stochastic optimal control problem such that a controller steers the covariance of a stochastic system to a target terminal value, while minimizing a state and control dependent cost in expectation.
Since the late `80s <cit.>, the infinite horizon case for linear time invariant systems has been thoroughly researched.
On the other hand, until recently, the finite horizon case has not been investigated <cit.>.
Our previous work <cit.> introduced state chance constraints into the OCS problem.
While it is possible to separately design mean steering (feedforward) and covariance steering (feedback) controllers, when state chance constraints exist, the state mean and state covariance constraints are coupled, and one needs to simultaneously design the mean and covariance steering controllers. 
In addition, in <cit.> we proposed an OCS controller that is computationally efficient and can deal with non-convex state chance constraints; the algorithm was then applied to vehicle path planning problems.

Note that the above mentioned OCS controllers cannot deal with input hard constraints, because they are affine functions of the state or the disturbance, which are both unbounded, and the control commands become unbounded as well.
Thus, they cannot satisfy input hard constraint specifications. 
Such a situation occurs in many real-world scenarios. 
For example, in aircraft <cit.> or spacecraft <cit.> control problems the control command values have physical restrictions such as minimum/maximum thrust.

The contribution of this article is the introduction of an OCS controller that can accommodate input hard constraints.
To the best of our knowledge, input hard constraints have not been discussed in the framework of OCS, while input soft constraints have been investigated by several prior works. 
These include the incorporation of a maximum value of the expectation of quadratic functions of the control inputs <cit.>, and the maximum probability of control constraint violation <cit.>. 


Our formulation of the input constraint is different than these approaches, in the sense that we directly impose hard constraints on the input.
Inspired by the approach in <cit.>, we propose to use a saturation function into our OCS controller <cit.> to impose input hard constraints. 

The remainder of this paper is organized as follows.
Section <ref> formulates the problem and introduces the input constrained OCS problem setup.

In Section <ref>, we introduce the newly developed input hard constrained OCS controller. 
Section <ref> validates the effectiveness of the proposed approach using numerical simulations. 
Finally, Section <ref> summarizes the results and discusses some possible future research directions. 

We use P ≻ (≽) 0 to denote that the matrix P is symmetric positive (semi)-definite.
Also, we use P > (≥) 0 to denote element-wise inequalities of the matrix P.
The trace of P is denoted with (P), and 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(P_0,…,P_N) denotes the block-diagonal matrix with block-diagonal matrices P_0,…, P_N.
v is the 2-norm of the vector v, and P is the induced 2-norm of the matrix P.
x ∼𝒩(μ, Σ) denotes a random variable x sampled from a Gaussian distribution with mean μ and (co)variance Σ.
[x] denotes the expected value, or the mean, of the random variable x. 
Finally, (·) is the error function.



§ PROBLEM STATEMENT

In this section we formulate the input-hard constrained OCS problem.



 §.§ Problem Formulation

We consider the following discrete-time stochastic linear time-varying system with additive noise

    x_k+1 = A_kx_k + B_ku_k + w_k,

where k denotes the time index, x∈^n_x is the state, u∈^n_u is the control input, and w∈^n_x denotes a random vector distributed according to a zero-mean white Gaussian noise with the properties,

    𝔼[w_k] = 0, 𝔼[w_k_1w_k_2^⊤] = 
    	 D_k_1D_k_2^⊤,     k_1 = k_2,
    
    	0,

In addition, A_k, B_k, and D_k are known system matrices having appropriate dimensions. 

Furthermore, we assume that

    𝔼[x_k_1w_k_2^⊤]=0,    0 ≤ k_1 ≤ k_2 ≤ N.

We also assume that the state and control inputs are subject to the constraints

    x_k ∈𝒳,   u_k ∈𝒰,

for all k ≥ 0, where 𝒳⊆^n_x and 𝒰⊆^n_u are convex sets containing the origin.
Throughout this work, we assume that the sets 𝒳 and 𝒰 are represented as an intersection of N_s and N_c linear inequality constraints, respectively, as follows

    𝒳   = ⋂_j = 0^N_s-1{x: α_x,j^⊤ x ≤β_x,j},
    𝒰   = ⋂_s = 0^N_c-1{u: α_u,s^⊤ u ≤β_u,s},

where α_x,j∈^n_x and α_u,s∈^n_u are constant vectors, and β_x,j∈ and β_u,s∈ are constant scalars.

Notice that, since the system noise in (<ref>) is possibly unbounded, the state is unbounded as well. 
Thus, we formulate the state constraints x_k ∈𝒳 probabilistically, as chance constraints,

    (x_k ∉𝒳) ≤ϵ,

where ϵ∈ (0, 1).
We keep the second set inclusion in (<ref>) since for the control inputs, hard constraints are preferable.
Using Boole's inequality, (<ref>) and (<ref>) can be satisfied if

    (α_x,j^⊤ x_k > β_x,j)    ≤ p_x,j,

for j=0,…,N_s-1, where p_x,j are pre-specified such that 

    ∑_j=0^N_s-1 p_x,j≤ϵ.


The initial state x_0 ∈^n_x is a random vector that is drawn from a normal distribution according to 

    x_0∼𝒩(μ_0,Σ_0),
 
where μ_0 ∈^n_x and Σ_0 ∈^n_x × n_x. We assume that Σ_0 ≽ 0.

Our objective is to design a control sequence {u_0,…,u_N-1} that steers the system state x_k to a target Gaussian distribution at time step N, that is,

    x_N ∼𝒩(μ_f, Σ_f),

where μ_f ∈^n_x, Σ_f ∈^n_x × n_x.
We assume that Σ_f ≻ 0 and wish to minimize the state and control expectation-dependent quadratic cost

    J(u_0,…,u_N-1) = 𝔼[∑_k=0^N-1x_k^⊤ Q x_k + u_k^⊤ R u_k],

where Q ≽ 0 and R ≻ 0.

In summary, we wish to solve the following finite-horizon optimal control problem

	
    J(u_0,…,u_N-1) = 𝔼[∑_k=0^N-1x_k^⊤ Q x_k + u_k^⊤ R u_k], 
       subject to
        x_k+1 = A_kx_k + B_ku_k + w_k, 
       (α_x,j^⊤ x_k > β_x,j) ≤ p_x,j, ∀ j ∈ [0,N_s-1], 
       α_u,s^⊤ u_k ≤β_u,s, ∀ s ∈ [0,N_c-1],
        x_N = x_f ∼𝒩(μ_f,Σ_f).




	System (<ref>) is assumed to be controllable in the absence of constraints and disturbances, that is, x_f is reachable from x_0 for any x_f ∈^n_x, provided that w_k = 0 for k = 0,…,N-1. This assumption implies that given any x_f ∈^n_x and x_0 ∈^n_x, there exists a sequence of control inputs {u_0,…,u_N-1} that steers x_0 to x_f in the absence of disturbances or any constraints.




§ PROPOSED APPROACH


This section introduces the proposed approach to solve Problem (<ref>).
Instead of (<ref>), we use the following equivalent form of the system dynamics

    X =  x_0 +  U+ W,

where X∈^(N+1) n_x, U∈^Nn_u, and W∈^Nn_x represents the concatenated state, input, and disturbance vector, respectively, e.g.,  X = [ x_0^⊤ x_1^⊤⋯ x_N^⊤ ]^⊤,
while the matrices ∈^(N+1)n_x× n_x, ∈^(N+1)n_x× Nn_u, and ∈^(N+1)n_x× Nn_x are defined accordingly <cit.>.
Note that 

	
    𝔼[x_0x_0^⊤]    = Σ_0 + μ_0μ_0^⊤,
    𝔼[x_0W^⊤]    = 0, 
    𝔼[WW^⊤]    =𝚋𝚕𝚔𝚍𝚒𝚊𝚐(D_0D_0^⊤,,…,D_N-1D_N-1^⊤).


We use the matrix E_k = [0_n_x,kn_x, I_n_x,0_n_x,(N-k)n_x]∈^n_x×(N+1)n_x and F_k = [0_n_u,kn_u, I_n_u,0_n_u,(N-1-k)n_u]∈^n_u× Nn_u such that x_k = E_k X and u_k = F_k U.



 §.§ Input-Constrained OCS Problem 


In this section, we discuss the OCS with input hard constraints. 

We start by introducing a relaxed form of Problem (<ref>) summarized in the following lemma.


	Given (<ref>), (<ref>), (<ref>), and (<ref>), Problem (<ref>) can be converted to the following form with a convex relaxation of the terminal covariance constraint.
	
		
    J(U) = 𝔼[X^⊤Q̅ X + U^⊤R̅ U ],
       subject to
       X =  x_0 +  U+ W,
       (α_x,j^⊤ E_k X > β_x,j) ≤ p_x,j,     j ∈ [0, N_s-1],  
       α_u,s^⊤ F_k U ≤β_u,s,  s ∈ [0, N_c-1], 
       μ_f = E_N𝔼[X],
       Σ_f ≽ E_N(𝔼[XX^⊤] - 𝔼[X]𝔼[X]^⊤) E_N^⊤,

	

	where Q̅ = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(Q,…,Q,0) ∈^(N+1)n_x×(N+1)n_x and R̅ = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(R,…,R) ∈^Nn_u× Nn_u.



	The procedure to relax Problem (<ref>) to (<ref>) is straightforward using the result of <cit.> and by relaxing the original terminal covariance constraint Σ_f = E_N(𝔼[XX^⊤] - 𝔼[X]𝔼[X]^⊤) E_N^⊤ to (<ref>).
	We would like to highlight that the Gaussian distribution can be fully defined by its first two moments. 
	This relaxation makes the terminal state less distributed than the target distribution. 
	In many real-world applications <cit.> it may be preferable to achieve smaller uncertainty in the terminal state. 


In order to solve Problem (<ref>), we propose a control law that is summarized in the following theorem, which is the main result of this paper.
This control policy is an extension of our previous controller <cit.> and is inspired from the approach in <cit.>. 


	The control law
	
    u_k = v_k + K_k z_k,

	where z_k is governed by the dynamics
	
		
    z_k+1   = A z_k + φ(w_k), 
    
    		z_0    = φ(ζ_0),   ζ_0 = x_0 - μ_0,
 
	
	where φ(·): ^d ↦^d is an element-wise symmetric saturation function with pre-specified saturation value of the ith entry of the input ζ_i^ max > 0 
	as 
	
    φ_i(ζ) = max(-ζ_i^ max, min(ζ_i,ζ_i^ max) ),

	converts Problem (<ref>) to the following convex programming problem
	
		
    J(V, K, Ω) =
       (Q̅[ I K ]Σ_XX[     I; K^⊤^⊤ ]) + (R̅KΣ_UU K^⊤) 
       + (μ_0 +  V)^⊤Q̅(μ_0 +  V) +  V^⊤R̅V 
       subject to 
       α_x,j^⊤ E_k(μ_0 +  V)- β_x,j +
       √(1-p_x,j/p_x,j)Σ_XX^1/2[ I K ]^⊤
    		E_k^⊤α_x,j≤ 0, 
       HF_kV + Ω^⊤σ≤ h,
       HF_kK[
    		 
    		] = Ω^⊤ S,
       Ω≥ 0,
       μ_f = E_N(μ_0 +  V),
       Σ_f ≽ E_N [ I K ]Σ_XX[     I; K^⊤^⊤ ] E_N^⊤,

	
	where Ω∈^2(N+1)n_x× N_c is a decision (slack) variable, 
	
    Σ_XX = 
       [  ;   ][              Σ_0    [ζ_0φ(ζ_0)^⊤];    [φ(ζ_0)ζ_0^⊤] [φ(ζ_0)φ(ζ_0)^⊤] ][ ^⊤   ;    ^⊤ ]
        + [  ;   ][       [WW^⊤]    [Wφ(W)^⊤];    [φ(W)W^⊤] [φ(W)φ(W)^⊤] ][ ^⊤   ;    ^⊤ ],

	
    Σ_UU = 
    		[φ(ζ_0)φ(ζ_0)^⊤]^⊤ + [φ(W)φ(W)^⊤]^⊤.

	Furthermore, 
	
		
    H    = [     α_u,0,         ⋯, α_u, N_c-1 ]^⊤∈^N_c× n_u, 
    
    		h    = [     β_u,0,          ⋯ β_u, N_c-1 ]^⊤∈^N_c,

	
	and
	
    V    = [   v_0^⊤       ⋯ v_N-1^⊤ ]^⊤,
    
    		K    = [ 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(K_0,K_1,…,K_N-1)              0_Nn_u× n_x ].

	
	In addition, S ∈^2(N+1)n_x×(N+1)n_x and σ∈^2(N+1)n_x are constant.
	Specifically, for i = 1,…, (N + 1)n_x
	
		
    S_2i-1   = e_2i-1^⊤,S_2i = -e_2i^⊤,
    σ_2i-1   = ζ_i^ max, σ_2i = ζ_i^ max.

	
	where S_i denotes the ith row of S, and e_i ∈^2(N+1)n_x is a unit vector with ith element 1.



	It follows from (<ref>) and (<ref>) that 
	
    Z = φ(ζ_0) + φ(W),

	where Z = [ z_0^⊤     … z_N^⊤ ]^⊤∈^(N+1)n_x.
	In addition, U is represented as U = V + K Z,
	and thus,
	
    U = V + K(φ(ζ_0) + φ(W)).

	Since the saturation function (<ref>) is an odd function and ζ_0 is zero-mean Gaussian distributed,
	
    [φ(ζ_0)] = 0,    [φ(W)] = 0,













	and hence,
	
    [U] = V,  Ũ = U - [U]
    		=K(φ(ζ_0) + φ(W)).

	Thus, it follows from (<ref>) that 
	
    X̅   = [X] = μ_0 +  V,
    X̃   = X -[X] 
    
    
    		   = [ I K ](
    		[  ;   ][    ζ_0; φ(ζ_0) ]
    		+ 
    		[  ;   ][    W; φ(W) ]).

	It also follows from (<ref>), (<ref>), (<ref>), and (<ref>) that 
	
    [X̃X̃^⊤]    = [ I K ]Σ_XX[     I; K^⊤^⊤ ],
    [ŨŨ^⊤]    =  KΣ_UU K^⊤.

	Following <cit.>, it can be shown that the cost function (<ref>) may be written as 
	
    J(V,Ũ)    = (Q̅[X̃X̃^⊤]) + X̅^⊤Q̅X̅
       + (R̅[ŨŨ^⊤]) + V^⊤R̅ V.

	Using (<ref>), (<ref>), and (<ref>), we can further convert (<ref>) to (<ref>).






	
	Next, we discuss the conversion of the chance constraint (<ref>) to a more amenable form to facilitate the numerical solution of the optimization problem (<ref>).
	First, notice that α_j^⊤ E_k X is a univariate random variable with mean α_x,j^⊤ E_k [X] and variance α_x,j^⊤ E_k [X̃X̃^⊤] E_k^⊤α_x,j.
	It follows from the Chebyshev-Cantelli inequality (Lemma <ref> in the Appendix) that 
	
    (α_x,j^⊤ E_k X ≤α_x,j^⊤ E_k [X] 
       +√(1-p_x,j/p_x,j)√(α_x,j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_x,j^⊤)) > 1-p_x,j.

	Therefore, the inequality (<ref>) is satisfied if
	
    α_x,j^⊤ E_k [X] +√(1-p_x,j/p_x,j)√(α_x,j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_x,j^⊤)≤β_x,j,

	is satisfied, which is, from (<ref>), (<ref>), and a similar discussion to <cit.>, equivalent to a second order cone constraint in terms of V and K (<ref>).
	

	In addition, using (<ref>) and (<ref>), the terminal state constraints (<ref>) and (<ref>) can be written as
	
		
    μ_f    = E_N(μ_0 +  V),
    Σ_f    ≽ E_N [ I K ]Σ_XX[     I; K^⊤^⊤ ] E_N^⊤.

	
	
	Using the Schur complement, (<ref>) can be further converted to 
	
    [                     Σ_f     E_N [ I K ]Σ_XX^1/2; Σ_XX^1/2[ I K ]^⊤ E_N^⊤                       I ]≽ 0.





	Note that Σ_XX≽ 0 (see Lemma <ref> in the Appendix).
	

	
	Finally, we rewrite the input hard constraint (<ref>) as follows.
	Using (<ref>), we first rewrite (<ref>) to the following equivalent form
	
    H F_k U ≤ h.

	Then, using (<ref>), this inequality is further converted to
	
    HF_k(V + K (𝒜φ(ζ_0) + 𝒟φ(W)) ) ≤ h,

	and thus, 
	
    HF_kV + HF_kK [ 𝒜 𝒟 ][ φ(ζ_0);   φ(W) ]≤ h.

	In addition, using (<ref>), the condition (<ref>) can be represented as 
	
    S [ φ(ζ_0);   φ(W) ]≤σ.

	It follows from the discussion in <cit.> that the constraint (<ref>) can be converted to 
	
	
		
    HF_kV + Ω^⊤σ≤ h,
       HF_kK[
    		 
    		] = Ω^⊤ S,
       Ω≥ 0.

	
	In summary, we have converted Problem (<ref>) to Problem (<ref>).


Note that the values of [ζ_0φ(ζ_0)^⊤], [φ(ζ_0)φ(ζ_0)^⊤], [Wφ(W)^⊤], and [φ(W)φ(W)^⊤] can be obtained using Monte Carlo or Lemma <ref> in the Appendix.


	When designing u_k in (<ref>), the value of z_k can be obtained from z_k-1 and φ(w_k-1).
	The noise w_k-1 is obtained from the system dynamics (<ref>), i.e., 
	
    w_k-1 = x_k - A_k-1x_k-1 - B_k-1u_k-1,

	and thus, φ(w_k-1) can be computed before computing the value of u_k at time step k. 




	Instead of (<ref>) or (<ref>), one can separately design the mean and covariance steering cost, i.e.,  	







	
    J(V,K) = (Q̅_v [ I K ]Σ_XX[ I K ]^⊤) + V^⊤R̅_mV 
       +(μ_0 +  V)^⊤Q̅_m(μ_0 +  V) + (R̅_vKΣ_UU K^⊤),

	where Q̅_m ∈ ^(N+1)n_x×(N+1)n_x, Q̅_v ∈ ^(N+1)n_x×(N+1)n_x, R̅_m ∈ ^Nn_u× Nn_u, and R̅_v ∈ ^Nn_u× Nn_u are block diagonal matrices, e.g., Q̅_m = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(Q_m,…, Q_m, 0). 



	Because the terminal covariance constraint (<ref>) can be formulated as a linear matrix inequality (LMI) constraint, in order to solve Problem (<ref>), we need to use an SDP solver such as MOSEK <cit.>.




§ NUMERICAL SIMULATION

In this section we validate the proposed OCS algorithm using a numerical example. 
We consider a path-planning problem for a vehicle under the following double integrator dynamics with x_k = [x,y,v_x,v_y]^⊤∈^4, u_k =[a_x,a_y]^⊤∈^2, w_k ∈^4 and

    A = [   1   0 Δ t   0;   0   1   0 Δ t;   0   0   1   0;   0   0   0   1 ], 
    	B = [ Δ t^2/2       0;       0 Δ t^2/2;     Δ t       0;       0     Δ t ],

and D = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(0.01, 0.01, 0.01, 0.01),
where Δ t is the time-step size set to Δ t = 0.2. 
The feasible state space is 

    0.2(x-1) ≤ y ≤ -0.2(x-1).

The mean and the covariance of the initial and target terminal distributions are set to

    μ_0 = [-10, 1, 0, 0], Σ_0 = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(0.05, 0.05, 0.01, 0.01),
       μ_f = [0, 0, 0, 0],  Σ_f = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(0.025, 0.025, 0.005, 0.005).

The cost function weights are chosen as 

    Q = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(0.5,4.0,0.05,0.05),  R = 𝚋𝚕𝚔𝚍𝚒𝚊𝚐(20,20).

The horizon is set to N=20, and the probability threshold to p_x,j = 0.05 for j = 0, 1.
Finally, we restrict the maximum acceleration at U_ max= 2.9 m/s^2 along each axis, i.e.,

    u_k _∞≤ U_ max,

for k= 0,…,N-1.
We also set the saturation function (<ref>) to be saturated when the input exceeds the 3σ values.
We employ YALMIP <cit.> along with MOSEK <cit.> to solve this problem.


We first show the results when the controller from <cit.> is used in Fig. <ref>.
The trajectories are depicted in Fig. *fig:Traj_F.
Red circles denote the initial and the target distributions, and blue ellipses represent the 3σ confidence region at each time step. 
We also illustrate randomly picked 100 trajectories with gray lines and observe that the state chance constraints are satisfied. 
Note, however, that the controller in <cit.> cannot deal with input hard constraints.
Figure *fig:ACC_F depicts the acceleration commands of the 100 sample trajectories with gray lines along with the acceleration limits ± 2.9 m/s^2 with red dashed lines. 
The cost for this scenario is 2,285.

Next, we show the results when the newly developed OCS controller in Theorem <ref> is used. 
Figure <ref> depicts the results. 
While, as shown in Fig. *fig:Traj_C, the state chance constraints are satisfied, the control commands depicted in Fig. *fig:ACC_C satisfy the input hard constraint (<ref>).
The cost for this constrained scenario is 2,301, which is, as expected, slightly larger than the cost for the unconstrained case owing to the additional input constraint (<ref>).







§ SUMMARY


In this work, we have addressed the problem of OCS under state chance constraints and input hard constraints.
Similarly to our previous works <cit.>, we solved this problem by converting the original problem into a convex programing problem.
The input hard constraints are formulated using saturation functions to limit the effect of possibly unbounded disturbance.
Numerical simulations show that the proposed algorithm successfully constructs control commands that satisfy the state and input constraints.
Future work includes OCS with measurement noise. 



IEEEtran




§ APPENDIX



	If a random variable z∈ is sampled from a Gaussian distribution 𝒩(0, σ_z^2) and the saturation function φ(·) is such that 
	
    φ(z) = max(-ζ, min(z, ζ)),

	where ζ > 0 is a constant, then 
	
    [φ(z)^2] =ζ^2 + (σ_z^2 - ζ^2)(ζ/√(2)σ_z) - 2ζσ_z/√(2π)e^-ζ^2/2σ_z^2.
       [zφ(z)] = σ_z^2(ζ/√(2)σ_z).



	The proof is straightforward from standard integral calculus. 



	The matrix Σ_XX defined as
	
    Σ_XX =
       [  ;   ][       [ζ_0ζ_0^⊤]    [ζ_0φ(ζ_0)^⊤];    [φ(ζ_0)ζ_0^⊤] [φ(ζ_0)φ(ζ_0)^⊤] ][ ^⊤   ;    ^⊤ ]
       + [  ;   ][       [WW^⊤]    [Wφ(W)^⊤];    [φ(W)W^⊤] [φ(W)φ(W)^⊤] ][ ^⊤   ;    ^⊤ ].

	is symmetric positive definite. 


	Consider vectors x ∈ ^(N+1)n_x and  y ∈ ^(N+1)n_x.
	Then, 
	
    [ x^⊤ y^⊤ ]Σ_XX[ x; y ]
       = [ x^⊤ y^⊤ ][       [ζ_0ζ_0^⊤]    [ζ_0φ(ζ_0)^⊤];    [φ(ζ_0)ζ_0^⊤] [φ(ζ_0)φ(ζ_0)^⊤] ][ ^⊤ x; ^⊤ y ]
        +[ x^⊤ y^⊤ ][       [WW^⊤]    [Wφ(W)^⊤];    [φ(W)W^⊤] [φ(W)φ(W)^⊤] ][ ^⊤ x; ^⊤ y ],
       =[ x^⊤ y^⊤ ][[    ζ_0; φ(ζ_0) ][    ζ_0^⊤ φ(ζ_0)^⊤ ]]
    	[ ^⊤ x; ^⊤ y ]
        + [ x^⊤ y^⊤ ][[    W; φ(W) ][    W^⊤ φ(W)^⊤ ]]
    	[ ^⊤ x; ^⊤ y ],


	
    =[(x^⊤ζ_0 + y^⊤φ(ζ_0))
    	(ζ_0^⊤^⊤ x + φ(ζ_0)^⊤^⊤ y)] 
        + [(x^⊤ W + y^⊤φ(W))
    	(W^⊤^⊤ x + φ(W)^⊤^⊤ y)],
       =[(x^⊤ζ_0 + y^⊤φ(ζ_0))^2] 
       + [(x^⊤ W + y^⊤φ(W))^2] ≥ 0.



[Chebyshev-Cantelli inequality <cit.>]
	Let z ∈ be a random variable with mean μ_z and variance Σ_z. Then, for any c≥ 0, the following inequality holds
	
    (z ≥μ_z + c) ≤Σ_z/Σ_z + c^2.
 

	
The following lemma describes the derivation of the chance constraints (<ref>) in detail. 


The chance constraint

    (α_x,j^⊤ E_k X > β_x,j) ≤ p_x,j,

can be satisfied if 

    α_x,j^⊤ E_k(μ_0 +  V)- β_x,j +
       √(1-p_x,j/p_x,j)Σ_XX^1/2[ I K ]^⊤
    		E_k^⊤α_x,j≤ 0




	It follows from Lemma <ref> that the following inequality on the random scalar variable α_j^⊤ E_k X holds
	
    (α_j^⊤ E_k X ≤α_j^⊤ E_k [X] + c) > 
        1 - α_j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_j^⊤/α_j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_j^⊤ + c^2,

	In addition, inequality (<ref>) can be equivalently rewritten as
	
    (α_x,j^⊤ E_k X ≤β_x,j) > 1 - p_x,j.

	We wish to compute c>0 that satisfies
	
    p_x,j = α_j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_j^⊤/α_j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_j^⊤ + c^2,

	and obtain 
	
    c = √(1-p_x,j/p_x,j)√(α_x,j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_x,j^⊤).

	Therefore, the following inequality holds
	
    (α_x,j^⊤ E_k X ≤α_x,j^⊤ E_k [X] 
       +√(1-p_x,j/p_x,j)√(α_x,j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_x,j^⊤)) > 1-p_x,j.

	By comparing (<ref>) and (<ref>), it follows that (<ref>) is satisfied if
	
    α_x,j^⊤ E_k [X] +√(1-p_x,j/p_x,j)√(α_x,j^⊤ E_k [X̃X̃^⊤]E_k^⊤α_x,j^⊤)≤β_x,j,

	which is equivalent to the second order cone constraint in terms of V and K (<ref>).


