P_l\!-\!\displaystyle\sum_{i:X_{li}=1}\!D_i^d(t)\!-\!\textbf{E}\!\{N_l(t)\!\}\!\geq\!\frac{\sigma_{N_l}}{\sqrt{2\epsilon}}.
\end{small}$$

Real-time Decomposition
-----------------------

In order to develop algorithm to regulate the SOC of each battery to satisfy (\[eq\_5\]), the shifted version of $B_i(t)$, $H_i(t)\!=\!B_i(t)\!-\!T_{max}\!-\!\eta^-\!J_ir_{i,max}$ is given, where $T_{max}$ is a constant which will be defined in Section . Like $B_i(t)$, $H_i(t)$ also satisfies $$\label{eq_10}
\begin{small}
H_{i}\!\left( t+1\!\right)\!=\!H_{i}\left( t\right)\!-\!\eta ^{-}\!\left(\!\left(\!\sum\limits_{j=1}^{J_{i}}r_{ij}\left( t\right)\!-\!U_{i}\left( t\right)\!\right)^{+}\!-\!D_{i}^{d}\left( t\right)\!\right)\!+\!\eta ^{+}R_{i}\left( t\right).
\end{small}$$   The technique for shifted queue is also utilized in [@Guo]. Then the state of system queues is described as $\overrightarrow{\theta(t)}=(Q_{11}(t),...,Q_{IJ_I}(t),H_1(t),...,H_I(t))$.

Define Lyapunov function as: $$\label{eq_11}
\begin{small}
\mathcal{L}(\overrightarrow{\theta(t)})=\frac{1}{2}\sum_{i=1}^I\sum_{j=1}^{J_i}Q_{ij}^2(t)+\frac{1}{2}\sum_{i=1}^IH_i^2(t).
\end{small}$$

The Lyapunov drift is given by $$\label{eq_12}
\begin{small}
\triangle(\overrightarrow{\theta(t)})=\textbf{E}\{\mathcal{L}(\overrightarrow{\theta(t+1)})-\mathcal{L}(\overrightarrow{\theta(t)})| \overrightarrow{\theta(t)}\}.
\end{small}$$

From (\[eq\_6\]) and (\[eq\_10\]), it is easy to see that $$\label{eq_13}
\begin{small}
 \frac{1}{2}(Q_{ij}^2(t+1)\!-\!Q_{ij}^2(t))\!\leq\!\alpha_{\max}\!-\!Q_{ij}(t)r_{ij}(t)
 \!-\!r_{ij}(t)\!\sum\nolimits_{k=1}^{K}\!E_k(t)W_{k,ij}(t)1_{\{Q_{ij}(t)=0\}},
 \end{small}$$ $$\label{eq_14}
\begin{small}
\frac{1}{2}(H_{i}(t+1)\!-\!H_{i}^2(t))\!\leq\!\beta_{\max}\!-\!H_{i}(t)\!\left(\eta^{-}\!\left(\!\left(\!\sum\limits_{j=1}^{J_{i}}r_{ij}(t)\!-\!U_{i}( t)\right)^{+}\!-\!D_{i}^{d}(t)\!\right)\!+\!\eta ^{+}R_{i}(t)\right),
\end{small}$$ where $\alpha_{\max}\!\!=\!\!\frac{1}{2}r_{max}^2\!\! +\!\! \frac{1}{2}E_{max}^2\!\! <\!\! \infty$ and $\beta _{\max }=\frac{1}{2}\max \left\{ \left( \eta ^{+}U_{\max }\right)^{2},\max_{i}\left\{ J_{i}^{2}r_{i,\max }^{2}\right\} \right\}\!\!<\!\!\infty$ are finite constants.

Substitute (\[eq\_13\]) and (\[eq\_14\]) into (\[eq\_12\]). We consider the following drift plus penalty expression as in [@Jin]. The Lyapunov drift is used to measure stability and the penalty item accounts for system performance. The tradeoff between system stability and optimal performance is achieved by the parameter $V$. $$\label{eq_15}
\begin{split}
&\triangle(\overrightarrow{\theta(t)})\!+\!V\textbf{E}\{c(t)\sum_{i=1}^{I}D_i^d(t)|\overrightarrow{\theta(t)}\}
\!\leq\!\delta_{\max}\!-\!\sum\limits_{i=1}^{I}\!\mathbf{E}\!\left\{\sum\limits_{j=1}^{J_{i}}r_{ij}\left( t\right)\!\left( Q_{ij}\left( t\right)
\!+\!H_{i}\left( t\right)\eta ^{-}\!\right)\!\cdot\!1_{\left\{ Q_{ij}\left(
t\right)>0\right\} }\left\vert \overrightarrow{\theta \left( 