

\subsection{Operation Procedure}
\begin{figure}[tbp]
    \centering
    \includegraphics[width=.75\columnwidth]{figs/flowchart.pdf}
    \caption{Real-time operation procedure}
    \label{fig:flowchart}
\end{figure}

This section introduces the real-time operation procedure for the DSO to manage the grid. PV capacity, EV power and energy capacity, HP capacity, water tank volume and temperature range are assumed known parameters to the DSO upon their installation. As shown in Fig. \ref{fig:flowchart}, before each operation step, the DSO collects power profiles from the customers via ICT infrastructures for the next time step, and runs a power flow program to evaluate the grid states. In case bus voltage limits or cable loading limits are violated, the DSO again collects data for model parameters and activates the coordination model to optimize the active and reactive power setpoints of various controllable DERs. It is assumed water tanks are used for domestic hot water heating and space heating. Their temperature is thus used as an important parameter to determine the urgentness of their power demand. Correspondingly, EVs' state of charge (SoC) and departure time are used to determine their urgentness. Finally, the optimized active and reactive power setpoints are re-fed into the power flow program for evaluation and grid state extraction considering that the proposed coordination model adopts relaxation for the power flow relations. Control actions are implemented afterwards.

\subsection{DER Coordination Model}
\begin{table}[tbp]
\begin{mdframed}
\input{nomenclature}
\printnomenclature
\end{mdframed}
\end{table}

This section presents the convex DER coordination model, which is implemented when grid limits are violated. Notations used are listed in the nomenclature. The model aims to restore grid limits by exploiting flexibility from various DERs, in which the decision variables include active and reactive power setpoints of EVs, HPs, and PVs connected to the grid. The objective function (\ref{obj}) minimizes the deviation from the DER users' original schedules collected from the customers, considering meanwhile losses in the grid. The cost terms are determined such that DER users with more urgent demand are assigned with higher values, which lowers the chance of their demand being curtailed. In this way, user satisfaction is more likely to be met. 

The constraints in the model are presented in the follows. The active power and reactive power injections into the grid from each household are defined with (\ref{pbalance}) and (\ref{qbalance}), respectively. The PVs' apparent power ratings are enforced in (\ref{pv1}), while (\ref{pv2}) reflects PV's power factor limits. Similar constraints for EVs are defined in (\ref{ev1}) and (\ref{ev2}). Fixed power factors for HPs are enforced in (\ref{hp}). Constraints (\ref{evcurtail}) and (\ref{pvcurtail}) respectively define active power curtailment for EVs and PVs. Constraints (\ref{hpcurtail1})-(\ref{hpcurtail3}) combined define the downward and upward regulation of active power for HPs.

The power flow relations are defined in (\ref{pflow})-(\ref{ilimits}), which adopts a second-order cone formulation \cite{Farivar2013}. Specifically, (\ref{pflow}) and (\ref{qflow}) define the active power and reactive power balance for each node, respectively. Voltage relation is enforced in (\ref{vrelations}), which is an exact formulation. The second-order cone relaxation is implemented in (\ref{socp}), which is relaxed from its equality form, defining apparent power flow for each cable. Constraints (\ref{vlimits}) and (\ref{ilimits}) define bus voltage limits and cable loading limits. Although there exist other convex relaxation models of the power flow relations such as semidefinite programming formulation \cite{Lavaei2012}, the second-order cone relaxation generally has a high accuracy and is exact under mild conditions for radial networks \cite{Farivar2013}, thus is widely used. 

\begin{subequations}
\begin{IEEEeqnarray}{rl} \label{obj}
&\underset{(p_j^{ev},p_j^{hp},p_j^{pv},q_j^{ev},q_j^{pv},\forall j \in \mathcal{H})}{\text{minimize}}\Big\{c^{loss}\sum_{\forall (i, j) \in \mathcal{E}} r_{i j} \ell_{i j} \nonumber\\&   +  \sum_{\forall j \in \mathcal{H}}(c_j^{ev} p_j^{ev\downarrow} + c_j^{pv} p_j^{pv\downarrow} + c_j^{hp\uparrow} p_j^{hp\uparrow}+ c_j^{hp\downarrow} p_j^{hp\downarrow} )\Big\} 
\end{IEEEeqnarray}
subject to
\begin{equation} \label{pbalance}
   p_j = p_j^{pv}-p_j^{load}-p_j^{ev}-p_j^{hp}  , \forall j \in \mathcal{H}
\end{equation}
\begin{equation} \label{qbalance}
   q_j = q_j^{pv}+q_j^{ev}-q_j^{hp}-q_j^{load} , \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{pv1}
    (p_j^{pv})^2 +(q_j^{pv})^2 \leq (S_j^{pv})^2, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{pv2}
    - p_j^{pv} \tan \Phi^{pv,max} \leq q_j^{pv} \leq p_j^{pv} \tan \Phi^{pv,max}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{ev1}
    (p_j^{ev})^2 +(q_j^{ev})^2 \leq (S_j^{ev})^2, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{ev2}
    - p_j^{ev} \tan \Phi^{ev,max} \leq q_j^{ev} \leq p_j^{ev} \tan \Phi^{ev,max}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{hp}
    q_j^{hp} = p_j^{hp} \tan \Phi^{hp}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{evcurtail}
    0 \leq p_j^{ev} \leq p_j^{ev,max}, p_j^{ev\downarrow} = p_j^{ev,max}-p_j^{ev}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{pvcurtail}
    0 \leq p_j^{pv} \leq p_j^{pv,fore}, p_j^{pv\downarrow} = p_j^{pv,fore}-p_j^{pv}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{hpcurtail1}
    0 \leq p_j^{hp} \leq p_j^{hp,max},  \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{hpcurtail2}
    p_j^{hp\downarrow} \geq \max(0, p_j^{hp,set}-p_j^{hp}), \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{hpcurtail3}
    p_j^{hp\uparrow} \geq \max(0, p_j^{hp}-p_j^{hp,set}), \forall j \in \mathcal{H}
\end{equation}
\begin{equation}\label{pflow}
p_{j}=\sum_{k: j \rightarrow k} P_{j k}-\sum_{i: i \rightarrow j}\left(P_{i j}-r_{i j} \ell_{i j}\right), \forall j \in \mathcal{N} 
\end{equation}
\begin{equation}\label{qflow}
q_{j}=\sum_{k: j \rightarrow k} Q_{j k}-\sum_{i: i \rightarrow j}\left(Q_{i j}-x_{i j} \ell_{i j}\right), \forall j\in \mathcal{N}
\end{equation}
\begin{IEEEeqnarray}{rl}\label{vrelations}
v_{j}=v_{i}-2\left(r_{i j} P_{i j}+x_{i j} Q_{i j}\right)+&\left(r_{i j}^{2}+x_{i j}^{2}\right) \ell_{i j},\\
&\forall(i, j) \in \mathcal{E} \nonumber
\end{IEEEeqnarray}
\begin{equation}\label{socp}
P_{i j}^{2}+Q_{i j}^{2} \leq \ell_{i j} v_i , \forall(i, j) \in \mathcal{E}
\end{equation}
\begin{equation}\label{vlimits}
   v_j^{min} \leq v_j \leq v_j^{max},  \forall j\in \mathcal{N}
\end{equation}
\begin{equation}\label{ilimits}
   0 \leq \ell_{i,j} \leq \ell_{i,j}^{max},
\forall(i, j) \in \mathcal{E}
\end{equation}
\end{subequations}

In summary, this model can restore grid limits by making full use of DER flexibility. It is noteworthy that this model requires only a one-step-ahead forecast, without the need for a multi-period setting and a multi-period-ahead forecast. The time-coupling DERs are coped with by deviating as little as possible from their original schedules while prioritizing DERs with the most urgent demand using dynamic cost terms.
\subsection{Case Description}
This section presents a case study to evaluate the proposed procedure and DER coordination model for real-time grid operations. The simulation is performed on a Dutch LV network, shown in Fig. \ref{fig:epse}. A 100\% penetration scenario for EVs, PVs, and HPs is assumed for 67 connected residential households seeing the energy transition.
\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/epse.png}
    \caption{Dutch LV network topology}
    \label{fig:epse}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/dataShow.pdf}
    \caption{Base load, PV, EV, HP profiles for a household for the selected day, starting from 8:00 with a 15-minute resolution}
    \label{fig:datashow}
\end{figure}

At the first step, a power flow program is implemented for a year with a 15-minute resolution. From which, a winter day with the most significant under-voltage and overloading issues is selected to evaluate the proposed real-time operation procedure and coordination model. The base load profiles are based on anonymous smart meter data from a Dutch DSO, while the PV generation profile is calculated from an annual solar irradiance dataset from the Royal Netherlands Meteorological Institute assuming a 4.5kW PV installation for each household. The HP consumption profiles are based on anonymous smart meter gas usage data from a Dutch DSO, assuming 70\% of gas is used for domestic heating \cite{Ruhnau2019}. Literature \cite{Brunner2016, Sichilalu2017} can be referred to for HP and water tank models. For EV charging, the Monto Carlo simulation is performed by assuming normally distributed EV arrival time and a fixed electricity demand of 7.5kWh per day. A 230V, 16A charger with a 90\% efficiency is assumed for EV charging. For illustration purpose, base load, PV, EV, and HP profiles for a random household on the selected winter day is presented in Fig. \ref{fig:datashow}, where the coincidence of uncontrolled EV charging and HP consumption during nighttime is likely to result in under-voltage and overloading issues in the grid.

\subsection{Simulation Flowchart}
\begin{figure}[tbp]
    \centering
    \includegraphics[width=.75\linewidth]{figs/simFlow.pdf}
    \caption{Simulation flowchart}
    \label{fig:simFlow}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/CDF_voltage.pdf}
    \caption{Bus voltage CDF}
    \label{fig:voltageCDF}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/CDF_loading.pdf}
    \caption{Cable loading CDF}
    \label{fig:loadingCDF}
\end{figure}

This section introduces the setting for the simulation on the selected day. The programs are implemented on a laptop with an Intel Core i5-9300H CPU running at 2.4GHz with 8GB RAM. The power flow program is implemented with Pandapower \cite{pandapower.2018}, while the optimization model is solved with Gurobi \cite{gurobi} with Python interface. The power flow program takes on average 1 second to solve, while the optimization model takes on average 4 seconds. It is noteworthy that a much higher resolution simulation could be performed under such short execution time, a 15-minute resolution simulation is still used considering the granularity of the available data.

The simulation process is described in Fig. \ref{fig:simFlow}, which starts with collecting various datasets as specified. At each simulation time step, a power flow program is implemented, followed by the DER coordination model in case any grid limit is violated. The optimized schedule from the coordination model is implemented, following which the EV charging and HP consumption profiles are updated for subsequent time slots to catch up on the deviation from the original schedules. For EVs, the charging process is extended in case the charging is curtailed. For HPs, the consumption is adjusted within their full capacity ranges to cancel out the deviation with the least possible delay. Finally, the EV state of charge (SoC) and water tank temperature are updated, which are used to calculate their corresponding cost terms in the DER coordination model and to evaluate user comfort.

\begin{table}[tbp]
    \centering
        \caption{Cost terms for the DER coordination model}
    \begin{tabular}{|c|c|c|} \hline
  \textbf{Cost term} & \textbf{Value} & \textbf{Note}\\ \hline
       $c^{loss}$  & 32\texteuro/MWh & Average day-ahead price \\ \hline
       $c^{pv}$  & 200\texteuro/MWh & Behind-the-meter\\\hline
       $c^{ev}$  & $\frac{c_0(1-SoC)}{(T_{max}-t+1)}$\texteuro/MWh & $c_0=1440, T_{max}=96$\\\hline
       $c^{hp\uparrow}$  & $\max(10,150\Delta T)$\texteuro/MWh & $\Delta T$ normalized \\\hline
       $c^{hp\downarrow}$  & $\max(10,-150\Delta T)$\texteuro/MWh & $\Delta T$ normalized \\\hline
    \end{tabular}
    \label{tab:cost}
\end{table}

For EVs, the cost terms are assumed to be functions of SoC and the remaining charging time. A unified departure time is assumed for the EVs, at the end of the simulation period. User satisfaction is determined by whether the EVs are fully charged before departure. For HPs, water tanks are assumed to supply hot water and space heating. A water tank of 1000L is assumed for households with no more than 5kW heating demand, while an additional 200L is assumed for each kW above 5kW \cite{Williams2012}. The water tanks are allowed for ${5}^{\circ} \mathrm{C}$ deviation up and down from the reference temperature, while not affecting user comfort \cite{Brunner2016, Sichilalu2017}. The cost terms are assumed as linear functions of water tank temperature deviation with a lower bound. User satisfaction is evaluated with respect to the water tank temperature limits. The used cost terms for this case study are listed in Table \ref{tab:cost}. It is noteworthy that in real-world applications, these cost terms can be derived through extensive simulation to guarantee the model's performance, or directly from end-users from a market perspective. With these model settings, the simulation results are presented in the following sections.

\subsection{Network Issues}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/comparison.pdf}
    \caption{Voltage magnitude, EV/HP power and EV reactive power generation for a time slot}
    \label{fig:compare}
\end{figure}

As described above, the major issues envisioned in this case study are under-voltage and overloading, which have been confirmed in Fig. \ref{fig:voltageCDF} and Fig. \ref{fig:loadingCDF} respectively. The reported grid operation states are both from the power flow programs, instead of the optimization model. Using the DER coordination model, both issues have been completely resolved, where the bus voltage stays higher than 0.9 pu and the cable loading percentage is no more than 70\%. Fig. \ref{fig:compare} presents the voltage magnitude, EV and HP power demand, and EV reactive power generation for a time slot before and after implementation of the DER coordination model, which suggests that the restoration of voltage magnitude is attributed to curtailed EV and HP power demand, and EV reactive power compensation.

\subsection{User Satisfaction}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/evsatisfaction.pdf}
    \caption{EV charging profiles for various households}
    \label{fig:evsatisfy}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/hpsatisfaction.pdf}
    \caption{Water tank temperature deviation for various households}
    \label{fig:hpsatisfy}
\end{figure}

In this study, two indicators are used to evaluate user satisfaction: whether EVs have been fully charged, and whether the temperature deviation of water tanks has exceeded the upper and lower limits. Fig. \ref{fig:evsatisfy} shows 3 different charging profiles, where the first EV's charging process is not interrupted, while the rest two have been curtailed during the charging process. However, all EVs have been fully charged at the end of the simulation. Fig. \ref{fig:hpsatisfy} presents the water tank temperature changes for various households. The temperature has been quickly restored to the normal state after a HP curtailment occurs. It is shown that all water tanks have not violated the temperature limits (${5}^{\circ} \mathrm{C}$ up and down). Combining these results, it is demonstrated that user satisfaction has been kept using this DER coordination model.


\subsection{SOCP Relaxation Accuracy}
\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/socp.pdf}
    \caption{Voltage magnitude and cable loading from both SOCP relaxation model and power flow program}
    \label{fig:SOCP}
\end{figure}

This section deals with the accuracy of the adopted SOCP relaxation of the power flow relations. Fig. \ref{fig:SOCP} shows the voltage magnitude and cable loading for one time slot from the DER coordination model and the power flow program. It is seen that the voltage has been precisely captured by the SOCP model, while the cable loading has been over-estimated, due to the relaxation of the apparent power flow relation. However, an important finding is that the SOCP relaxation model ensures the feasibility of the optimized active and reactive power setpoints, which implies that when these control actions from the DER coordination model are implemented in the grid, it would not result in loading or voltage issues that are not captured in the SOCP model.



\section{Introduction}
\input{1.intro}
\section{DER Coordination model}\label{sec:model}
\input{2.mathmodel}
\section{Case study}\label{sec:case}
\input{3.case}
\section{Conclusions}\label{sec:conclusion}
\input{4.conclusions}













