\section{Introduction}
	{W}{ith} increasing concerns about energy shortages and global warming, integrating more renewable generations into the power system has become a worldwide trend. However, uncertain and intermittent renewable generations have brought new challenges to power systems' secure and reliable operations. To deal with these challenges, system operators need to properly procure reserve from dispatch-able resources to defend against possible contingencies and load/renewable generation fluctuations. Since energy and reserve are tightly coupled, all the independent system operators (ISOs) in the U.S. run joint optimizations \cite{review1}, by and large with the following model:
	\begin{align}
	\label{obj1}
	& \mbox{(I)}: \quad\underset{\{g,r_U,r_D\}}{\rm minimize} \quad C_{g}^T g + C_{U}^T r_U + C_{D}^T r_D\\
	&\mbox{subject to}\notag\\
	\label{traditional balance}
	&\lambda : \mathds{1}^T g =  \mathds{1}^T d,\\
	\label{traditional pf}
	&\mu : S(g-d)\leq f,\\
	\label{reserve requirement}
	&(\gamma^U,\gamma^D):\mathds{1}^T r_U =  R^U,\mathds{1}^T r_D =  R^D,\\
	\label{physical limit}
	&g + r_U \leq \overline{G}, \underline{G} + r_D \leq g,0 \leq r_U \leq \overline{r_U}, 0 \leq r_D \leq \overline{r_D}.
	\end{align}
	The objective function (\ref{obj1}) aims to minimize the bid-in cost of energy and reserve with decision variables $(g,r_U,r_D)$, which denote vectors of energy, upward reserve, and downward reserve procured from generators. Coefficients $(C_{g},C_{U},C_{D})$ denote bid-in price vectors of $(g,r_U,r_D)$. Constraints (\ref{traditional balance})-(\ref{physical limit}) represent, respectively, energy balancing constraints, line capacity limits, reserve requirement constraints, and generator capacity and ramping rate limits. Reserve clearing prices will be obtained from the Lagrangian multipliers $(\gamma^U,\gamma^D)$, representing the increase of total bid-in cost when there is one additional unit of upward/downward reserve requirement.

	Regarding the traditional co-optimization model (I), there are several key issues that remain unclear. The first is about reserve requirements $R^U$ and $R^D$. These parameters are artificial and empirical, sometimes specified as the capacity of the biggest generator as in PJM\cite{PJM} or a certain proportion of system loads as in CAISO\cite{CAISO}. However, they will significantly affect the clearing results and prices of reserve and even the energy market. The second is on the locations of reserve resources. The traditional model only considers network limits in the base case. When there are contingencies or load/renewable generation fluctuations, the deliverability of reserve may still be limited by line capacities. A common solution is to partition the entire system into different zones and specify zonal reserve requirements as in ISO-NE\cite{ISO-NE}, which is again empirical. The third is about the objective function (\ref{obj1}), which is to minimize the base-case bid-in cost but does not consider possible re-dispatch costs under contingencies and fluctuations.

	\indent In light of these problems, there are also many research efforts from academia. Some researches utilize the forecast methods, such as the density forecasts which consider the probability distributions of future observations\cite{MatosBessa1}, or the scenario forecasts which consider several typical future scenarios\cite{scen1}, to define reserve requirements in systems with high renewable penetrations, see \cite{requirementReview} for a survey. In \cite{reservezone}, a statistical clustering algorithm is proposed to enable dynamic reserve zone partition. In \cite{readjust1} and \cite{readjust2}, the energy balancing and network constraints considering events and re-dispatches in contingency scenarios are analyzed, and the reserve pricing approach at the bus level is proposed. In \cite{UMP}, the uncertainty marginal prices are defined to price both uncertainty sources and reserve resources at the bus level.

	\indent In this paper, a scenario-oriented energy-reserve co-optimization model is developed. These scenarios will correspond to possible contingencies or load/renewable generation fluctuations, or their combinations. Reserve procured from each generator is equal to the maximum range of its generation re-dispatch among all scenarios, thereby eliminating the need to specify the parameters of reserve requirements artificially. Transmission capacity limits in all scenarios are also included in the proposed model, therefore the artificial partition of reserve zones is no longer needed either. Thereupon, we derive marginal prices of generations, loads and reserve. Under certain assumptions, energy and reserve marginal prices will be locational uniform prices. We have also established the property of revenue adequacy for SOs: revenue from load payments, credits to generators including energy, reserve, and re-dispatches, and congestion rent will reach their balance in the base case as well as in all scenarios considered.

	\section{Energy-Reserve Co-optimization Model}
     To properly address the problems of the traditional model (\uppercase\expandafter{\romannumeral1}), a scenario-oriented energy-reserve joint optimization model is proposed. Consider the current market structure, the proposed model can be regarded as an intermediate look-ahead optimization stage between the day-ahead (DA) scheduling and the real-time security constrained economic dispatch (RT SCED), like the ancillary service optimizer (ASO) in PJM. As Fig. \ref{PJM RT SCED} shows, ASO is one of many look-ahead stages to procure reserve resources with different flexibility levels. In practice, ASO is only financial binding for reserve\cite{PJM}. In this paper, however, we consider the proposed co-optimization model to be financial binding for both energy and reserve.

	\begin{figure}[H]
		\centering
		\vspace{-0.35cm}
		\includegraphics[width=3.4in]{PJM_RT_SCED_2.pdf}
    \setlength{\abovecaptionskip}{-3.0pt}
    \setlength{\belowcaptionskip}{-2.5pt}
    \captionsetup{font=footnotesize,singlelinecheck=false}
		\caption{PJM real-time operations including ASO, IT SCED, RT SCED and LPC\cite{PJM}}
		\label{PJM RT SCED}
    \vspace{-0.15cm}
	\end{figure}

	The proposed model makes the following assumptions:

	1) A standard DCOPF model with linear cost functions for energy and reserve is adopted, which is consistent with the current electricity market design.

	2) A single-period problem is considered for simplicity\footnote{A multi-period setting will bring in the coupling between ramping and reserve, which is a highly complicated issue.}.

	3) For simplicity, we don't consider the DA market. If the DA market clearing results need to be considered, the quantities in the proposed model can be seen as the deviations from the quantities in the DA scheduling, and all the qualitative analyses in this paper will still hold.

	4) Uncertainties from line outages and load fluctuations are considered. Generator outages are not considered\footnote{Generator outages will bring in the complicated issue of non-uniform pricing. Essentially, energy and reserve procured from generators with different outage probabilities are no longer homogeneous goods. We will skip this tricky problem now and leave it to our future work.}.

	5) Renewable generations are modeled as negative loads.

	6) The SO's predictions of possible scenarios as well as their occurrence probabilities are perfect.

	Based on assumptions (1)-(6), the proposed model is:
	\begin{equation}
	\setlength\abovedisplayskip{1pt
    \setlength\belowdisplayskip{-2.5pt}
    \begin{split}
	&\!F(g,r_U,r_D,\delta g^U_k,\delta g^D_k,\delta d_k)=\\
	\label{obj2}
	&\!C_{g}^T\!g\!+\!C_{U}^T\!r_U\!+\!C_{D}^T r_D\!+\!\sum^{k\in \mathcal{K}} \!\epsilon_k (\overline{C}^T_k\!\delta g^U_k \!-\!\underline{C}_k^T\!\delta g^D_k \!+\!C^T_{L}\!\delta d_k)\!,
    \end{split}
    \end{equation}
    \begin{align}
	&\mbox{(II)}: \quad \underset{\{g,r_U,r_D,\delta g^U_k,\delta g^D_k,\delta d_k\}}{\rm minimize} F(\cdot), \notag\\
	&\mbox{subject to}\notag\\
	\label{base balance and pf}
	&(\lambda,\mu):\mathds{1}^T g =  \mathds{1}^T d, S(g-d)\leq f,\\
	\label{base physical limit}
	&g + r_U \leq \overline{G}, \underline{G} + r_D \leq g,0 \leq r_U \leq \overline{r_U}, 0 \leq r_D \leq \overline{r_D},\\
	&\mbox{for all $k\in \mathcal{K}$:}
	\notag\\
	\label{cntg balance}
	&\lambda_k:\mathds{1}^T(g+\delta g^U_k-\delta g^D_k)=\mathds{1}^T (d+\pi_k -\delta d_k),\\
	\label{cntg pf}
	&\mu_k:S_k\bigl((g+\delta g^U_k-\delta g^D_k)-(d+\pi_k -\delta d_k)\bigr) \leq f_k,\\
	\label{dg1fanwei}
	&(\underline{\alpha_k},\overline{\alpha_k}): 0 \leq \delta g^U_k \leq r_U,\\
	\label{dg2fanwei}
	&(\underline{\beta_k},\overline{\beta_k}): 0 \leq \delta g^D_k \leq r_D,\\
	\label{dk1fanwei}
	&(\underline{\tau_k},\overline{\tau_k}): 0 \leq \delta d_k \leq d+\pi_k,
	\end{align}
    where $\mathcal{K}$ denotes the set of all non-base scenarios, coefficient $\epsilon_k$ denotes the occurrence probability of scenario $k$, decision variables $(\delta g^U_k,\delta g^D_k,\delta d_k)$ denote vectors of generation upward re-dispatches, downward re-dispatches, and load shedding in scenario $k$. Coefficients $(\overline{C_k},\underline{C_k})$ represent vectors of generation upward and downward re-dispatch prices\footnote{In some ISOs in the U.S., generation re-dispatches will be settled at the real-time LMPs plus some premiums or adders e.g. ERCOT\cite{ERCOT}.} in scenario $k$. If one generator's downward reserve is deployed, it will pay back to the SO, so there is a negative sign before $(\underline{C}_k^T\delta g^D_k)$ in (\ref{obj2}). Coefficient $C_L$ denotes the vector of load shedding prices. Coefficient $\pi_k$ denotes the vector of load fluctuations in scenario $k$. The objective function (\ref{obj2}) aims to minimize the expected system total cost, including the base-case bid-in cost and the expectation of re-dispatch costs in all scenarios. Constraints (\ref{base balance and pf})-(\ref{base physical limit}) denote the base-case limits as constraints (\ref{traditional balance})-(\ref{traditional pf}) and (\ref{physical limit}) in the traditional model (I). Constraints (\ref{cntg balance})-(\ref{cntg pf}) denote the energy balancing constraints and the network constraints in all non-base scenarios where $S_k$ and $f_k$ denote the shift factor matrix and the line capacity vector in scenario $k$. In non-base scenarios, the power flows on transmission lines are allowed to exceed line capacities for short duration, and the exceeding ratings are reflected in the setting of $f_k$. Constraints (\ref{dg1fanwei})-(\ref{dg2fanwei}) denote that generation re-dispatches in all scenarios shall not surpass the procured reserve, indicating that the procured reserve will be modeled as the maximum range of generation re-dispatch among all scenarios. Therefore reserve will be optimally procured wherever necessary to minimize the expected total cost. Constraint (\ref{dk1fanwei}) denotes that load shedding in each scenario must be non-negative and cannot exceed actual load power in that scenario.

	With such formulations, the solution to the proposed model will procure reserve at the best locations considering all kinds of costs. On the contrary, the solution to the traditional model will be either sub-optimal or infeasible in order to minimize the expected total cost subject to energy balancing and line capacity constraints in both the base case and all scenarios, even with reserve requirements $R^U$ and $R^D$ being the sum of the optimal procurement $r_U^*$ and $r_D^*$ in model (II).

	Moreover, marginal prices derived from the proposed scenario-oriented optimization model (II) have many attractive properties, as in the next section.

	\section{Pricing and Settlement}
	Based on the proposed model, in this section we present the pricing approach for generations, loads, and reserve, as well as the associated market settlement process.
	\subsection{LMP calculations for generations, loads and reserve}
	The proposed pricing approach for energy and reserve is based on their marginal contributions to the expected system total cost presented in (\ref{obj2}). Namely, consider any generator $j$, we first fix $g(j),r_U(j),r_D(j)$ at their optimal values $g(j)^*,r_U(j)^*,r_D(j)^*$ and consider them as parameters instead of decision variables. Such a modified optimization model is referred to as model (III). Compared with the original model (II), in model (III) the terms that are only related to generator $j$, specifically the bid-in cost of generator $j$ in the objective function (\ref{obj2}) and the $j^{th}$ row of all constraints in (\ref{base physical limit}), will be removed. Next, we evaluate the sensitivity of the optimal objective function of model (III), which represents the expected cost of all other market participants except for generator $j$, with respect to the deviations in parameters $g(j),r_U(j)$ and $r_D(j)$. According to the envelop theorem, we have:
	\begin{equation}
	\setlength\abovedisplayskip{-0.5pt
    \setlength\belowdisplayskip{2pt}
	\frac{\partial C^*_{\uppercase\expandafter{\romannumeral3}}}{\partial g(j)}\!=\!\frac{\partial \mathcal{L}_{\uppercase\expandafter{\romannumeral3}}}{\partial g(j)}\!=\!S(:,m_j)^T\mu\!-\!\lambda\!+\!\sum^{k\in \mathcal{K}} (S_k (:,m_j)^T \mu_k\!-\!\lambda_k),
	\end{equation}
	where $C^*_{\uppercase\expandafter{\romannumeral3}}$ denotes the optimal objective function of model (\uppercase\expandafter{\romannumeral3}) and $\mathcal{L}_{\uppercase\expandafter{\romannumeral3}}$ denotes the Lagrangian function of model (\uppercase\expandafter{\romannumeral3}). Accordingly the marginal energy price of any generator is:
	\begin{align}
	\label{RGMP2}
	\eta^g(j)&\!=\!- \frac{\partial \mathcal{L}_{\uppercase\expandafter{\romannumeral3}}}{\partial g(j)}\!=\!\lambda\!-\!S(:,m_j)^T\!\mu\!+\!\sum^{k\in \mathcal{K}}\!(\lambda_k\!-\!S_k(:,m_j)^T\!\mu_k)\notag\\
	&=\omega_0(j)+\sum^{k\in \mathcal{K}}\omega_k(j),
	\end{align}
	where $\eta^g(j)$ denotes the marginal energy price of generator $j$, with the term $m_j$ denoting bus $m$ where generator $j$ is located. The term $(\lambda-S^T\mu)$, denoted by $\omega_0$ for brevity hereafter, denotes the base-case component of the energy prices of generators. The term $(\lambda_k-S_k^T\mu_k)$, denoted by $\omega_k$ for brevity hereafter, denote the non-base component of the energy prices of generators in scenario $k$. We can see that each component has an energy part and a congestion part, as in the standard LMP formulation. Similarly, the marginal energy prices of loads are calculated as
	\begin{align}
	\label{L/d}
	\eta^d&\!=\!\frac{\partial C^*_{\uppercase\expandafter{\romannumeral3}}}{\partial d}\!=\!\frac{\partial \mathcal{L}_{\uppercase\expandafter{\romannumeral3}}}{\partial d}\!=\!\lambda\!-\!S^T\mu\!+\!\sum^{k\in \mathcal{K}}(\lambda_k\!-\!S_k^T\mu_k)\!-\!\sum^{k\in \mathcal{K}}\overline{\tau}_k
	\notag\\&\!=\!\omega_0+\sum^{k\in \mathcal{K}}\omega_k-\sum^{k\in \mathcal{K}} \overline{\tau}_k,
	\end{align}
	which is consistent with the marginal energy prices of generators except for the last term $(-\sum \overline{\tau}_k)$, which are multipliers associated with the upper bound of load shedding in (\ref{dk1fanwei}).

	According to the envelop theorem, there is
	\begin{equation}
	\label{RUMP}
	\eta^U(j)=-\frac{\partial C^*_{\uppercase\expandafter{\romannumeral3}}}{\partial r_U(j)}=\frac{\partial\mathcal{L}_{\uppercase\expandafter{\romannumeral3}}}{\partial r_U(j)}=\sum^{k\in \mathcal{K}} \overline{\alpha_k}(j),
	\vspace{-0.25cm}
	\end{equation}
	where $\eta^U(j)$ denotes the upward reserve marginal price of generator $j$. Similarly the downward marginal reserve price is:
	\begin{equation}
	\label{RDMP}
	\setlength\abovedisplayskip{-0.5pt
    \setlength\belowdisplayskip{1.5pt}
	\eta^{D}(j)=\sum^{k\in \mathcal{K}}\overline{\beta_k}(j).
	\end{equation}
	With such formulations, the connections between reserve and generation re-dispatches can be established: if the upward generation re-dispatch $\delta g^U_k(j)\!=\!r_U(j)$ in scenario k, then its corresponding multiplier $\overline{\alpha_k}(j)$ will be positive and contribute to the upward reserve marginal price $\eta^U(j)$.

	\indent Although the energy prices of generators and loads are defined separately, and the reserve marginal prices in (\ref{RUMP})-(\ref{RDMP}) are defined at the resource level, next we establish the uniform pricing property with some additional assumptions:

    7) Assume that $(-\sum\!\overline{\tau}_k)$ is zero in (\ref{L/d})\footnote{It is rare to shed a load to its total capacity. If a load is completely shed and another load at the same bus is not, this indicates that they have different reliability requirements, which will again bring in the non-uniform pricing issue.}.

	8) Assume that generators at the same bus will have the same $(\overline{C}_k,\underline{C}_k)$ in each scenario\footnote{The upward and downward re-dispatch prices $(\overline{C}_k,\underline{C}_k)$ in scenario $k$ can be set as the forecast RT-LMPs in scenario $k$ plus a fixed adder. With perfect predictions, the forecast RT-LMPs will be equal to the actual RT-LMPs, therefore the upward/downward re-dispatch prices $(\overline{C}_k,\underline{C}_k)$ can be used to settle the generation upward/downward redispatches.}.

	\begin{theorem}[Uniform Pricing]
	    Consider any two generators $i,j$ and any load $l$ at the same bus. Under assumptions (1)-(7), there is $\eta^g(i)=\eta^g(j)=\eta^d(l)$ for energy.

	   	Moreover, under assumptions (1)-(8) and assume that $r_U(i),r_D(i),r_U(j),r_D(j)\!>\!0$, there are $\eta^U\!(i)\!=\!\eta^U\!(j)$ for upward reserve and $\eta^D\!(i)\!=\!\eta^D\!(j)$ for downward reserve.
	\end{theorem}

	 The proof is presented in the Appendix A. Due to the page limit, we leave the proofs for this Theorem and the following Theorem 2 to the online version of our manuscript \cite{shi2020scenariooriented}.

	 The settlement process based on the proposed pricing approach has some attractive properties, as in the next subsection.

	\subsection{Market settlement process}
	Next the settlement process will be presented. The process can be separated into two stages. In the ex-ante stage, we don't know which prediction is true; and in the ex-post stage, the settlement depends on which scenario is actually realized.
	\subsubsection{ex-ante stage}
	 In this stage, generations and reserve will be financial binding and settled, also the basic load capacities $d$ and load fluctuations $\pi_k$ in all scenarios will be settled. Therefore, the ex-ante stage includes the following payments:
	\begin{itemize}
		\itemsep=2.2pt
		\item contribution of base-case prices to generator energy credit:
		\begin{equation}
		\label{base-case prices to generator energy credit}
		\setlength\abovedisplayskip{1pt}
        \setlength\belowdisplayskip{-0.5pt}
		\Gamma^g_0=\omega_0^Tg;
		\end{equation}

		\item contributions of non-base scenario prices to generator energy credit:
		\begin{equation}
		\label{non-base scenario prices to generator energy credit}
		\setlength\abovedisplayskip{1pt}
        \setlength\belowdisplayskip{-0.5pt}
		\sum\Gamma^g_k=\sum^{k\in \mathcal{K}}\omega_k^Tg;
		\end{equation}

		\item contributions of base-case price to load energy payment:
		\begin{equation}
		\label{base-case prices to load payment}
		\setlength\abovedisplayskip{-0.5pt}
        \setlength\belowdisplayskip{-0.5pt}
		\Gamma^d_0 = (\omega_0)^Td;
		\end{equation}

		\item contributions of non-base scenario prices to load energy payment:
		\begin{equation}
		\label{non-base scenario prices to load payment in all scenarios}
		\setlength\abovedisplayskip{-2pt}
		\setlength\abovedisplayskip{-1pt}
		\sum^{k\in \mathcal{K}}\Gamma^d_k = \sum^{k\in \mathcal{K}}(\omega_k)^T d;
		\end{equation}

	\item non-base load fluctuation payment in all scenarios:
		\begin{equation}
		\label{fluctuation payment}
		\setlength\abovedisplayskip{1.6pt}
		\setlength\abovedisplayskip{-1pt}
			\sum^{k\in \mathcal{K}}\Pi_k = \sum^{k\in \mathcal{K}}(\omega_k)^T \pi_k;
		\end{equation}

	\item upward and downward reserve credit:
		\begin{equation}
		\label{upward reserve payments}
	    \setlength\abovedisplayskip{0.5pt}

		\Gamma^U=(\eta^U)^T r_U=\sum^{k\in \mathcal{K}}\overline{\alpha_k}^T r_U=\sum^{k\in \mathcal{K}}\Gamma^U_k,
		\end{equation}
		\begin{equation}
		\label{downward reserve payments}
		\setlength\belowdisplayskip{-0pt}
		\Gamma^D=(\eta^D)^T r_D=\sum^{k\in \mathcal{K}}\overline{\beta_k}^T r_D=\sum^{k\in \mathcal{K}}\Gamma^D_k.
		\end{equation}
	\end{itemize}

	For the fluctuation payment, note that one load should pay for its fluctuations in all scenarios because the SO will procure reserve to deal with all possible load fluctuations accordingly, therefore the load fluctuation payments should not only depend on the realized scenario in real-time, but should also rely on other scenarios considered.
    \subsubsection{ex-post stage}
	In this stage, one of the non-base scenarios is actually realized and denoted as scenario $k$. With perfect predictions, for each generator $j$, the deviation of its real-time generation level from its base-case energy procurement $g(j)^*$ will be $\delta g^U_k(j)$ or $\delta g^D_k(j)$. Under the assumption (7) and the corresponding footnote (5), upward and downward generation re-dispatches will be settled with $\overline{C_k}$ and $\underline{C_k}$, and load shedding will be settled with the shedding prices $C_L$. Therefore, the ex-post stage includes the following payments:
	\begin{itemize}
	    \itemsep=1.0pt
		\item upward redispatch compensation:
		\begin{equation}
		\label{upward redispatch payments}
		\setlength\abovedisplayskip{-0.5pt}
		\setlength\belowdisplayskip{-1pt}
		\Phi^U_k\!=\!\overline{C}^T_k\delta g^U_k;
		\end{equation}
		\item downward redispatch pay-back:
		\begin{equation}
		\label{downward redispatch payments}
		\setlength\abovedisplayskip{-0.5pt}
		\setlength\belowdisplayskip{-1pt}
		\Phi^D_k\!=\!\underline{C}^T_k\delta g^D_k;
		\end{equation}
		\item load shedding credit:
		\begin{equation}
		\label{load shedding credits}
		\setlength\abovedisplayskip{-0.5pt}
		\setlength\belowdisplayskip{-1pt}
		\Phi^d_k\!=\!C_L ^T \delta d_k.
		\end{equation}
	\end{itemize}
	Note that the net revenue of the SO should be equal to the congestion rent. For the settlement process, the following theorem regarding SO's revenue adequacy is established:
	\begin{theorem}[Revenue Adequacy]
    Under assumptions (1)-(8), the net revenue of the SO is always non-negative.

	In particular, for the base case, load energy payment (\ref{base-case prices to load payment}) is equal to the sum of generator energy credit (\ref{base-case prices to generator energy credit}) and congestion rent:
	\begin{equation}
	\setlength\abovedisplayskip{-0.5pt}
	\setlength\belowdisplayskip{-0.5pt}
	\Gamma^d_0=\Gamma^g_0+\Delta_0=\Gamma^g_0+f^T\mu.
	\end{equation}

	For each non-base scenario k, its contribution to load payment is equal to the sum of its contribution to generator credit, load shedding credit, and congestion rent:
	\begin{equation}
	\setlength\abovedisplayskip{0.5pt}
	\setlength\belowdisplayskip{-0.5pt}
	\Gamma^d_k +\Pi_k=\Gamma^g_k+\Gamma^U_k+\Gamma^D_k+\epsilon_k\Phi^U_k -\epsilon_k\Phi^D_k+\epsilon_k\Phi^d_k+\Delta_k,k \in \mathcal{K},
	\end{equation}
	where the left-hand side is the contribution of scenario k to load payment, including energy payment (\ref{non-base scenario prices to load payment in all scenarios}) and load fluctuation payment (\ref{fluctuation payment}). The first five terms on the right-hand side represent its contribution to generator credit, including energy credit (\ref{non-base scenario prices to generator energy credit}), upward and downward reserve credit (\ref{upward reserve payments}-\ref{downward reserve payments}), upward and downward re-dispatch payment (\ref{upward redispatch payments}-\ref{downward redispatch payments}). The sixth term is the expected load shedding credit (\ref{load shedding credits}), and the last term is the congestion rent in that scenario $\Delta_k=f_k^T \mu_k$.
	\end{theorem}

    Please refer to the Appendix B in \cite{shi2020scenariooriented} for the proof. With this Theorem, payments from loads, payments to generators and congestion rent will reach their balance in the basecase as well as in all scenarios. Also, the reserve costs and the re-dispatch costs will be allocated in a scenario-oriented way.

	\section{Case Study}
    The case studies are performed both on a 2-bus system and on the modified IEEE 118-bus system.
	\subsection{Two-bus System}
	 \begin{figure}[H]
		\centering
		\captionsetup{font=footnotesize,singlelinecheck=false}
	    \vspace{-0.4cm}

		\includegraphics[width=2in]{one_line_diag.pdf}

		\caption{the One-Line Diagram for the 2-Bus System}
		\label{one line diag}
	    \vspace{-0.3cm}
	\end{figure}
    We first consider a 2-bus system with its one-line diagram presented in Fig. \ref{one line diag}. The generator bids are presented in Table \ref{Generators' Offers}. The outage probability of the two-parallel-line branch is 10\%, meaning that the system will lose one of these two lines if the outage happens. The line capacity exceeding rates are set to be for all scenarios. The basic loads are $(6,15,4)$MW, with two possible fluctuation situations: situation I $(+2,+6,-1)$MW with probability 20\%, and situation II $(+3,+2,-3)$MW with probability 20\%. Based on the line outage and load fluctuation information, we present all possible scenarios in Table \ref{Non-base}.

    	\begin{table}[H]
	    \footnotesize
		\centering
		\vspace{-0.2cm}
		\setlength{\abovecaptionskip}{-0.5pt}
        \setlength{\belowcaptionskip}{0pt}
        \captionsetup{font=footnotesize,textfont=sc}
		\caption{Generators' Offer Data for the 2-Bus System}

		\begin{tabular}{cccc}
			\hline
			Generator &$\overline{G}/\underline{G}$& $\overline{r_U}/\overline{r_D}$& $C_{g}/C_{U}/C_{D}$ \\
			\hline
		 G1  & $16/0$ & $4/4$ & $8/2/2$\\
		 G2  &$18/0$  &  $4/4$ & $15/2/2$  \\
		 G3 & $12/0$  &  $4/4$ &  $20/2.5/2.5$   \\
			\hline

			\end{tabular}
		\label{Generators' Offers}
		\vspace{-0.2cm}
		\end{table}

	\begin{table}[H]
	\centering
	\vspace{-0.3cm}
	\footnotesize
    \setlength{\abovecaptionskip}{-0.5pt}
    \setlength{\belowcaptionskip}{0pt}
        \captionsetup{font=footnotesize,textfont=sc}
			\caption{Non-Base Scenario Data for the 2-Bus System}
			\label{Non-base}
			\begin{tabular}{ccccc}
				\hline
				NO.  & Outage & Load Situation & Probability & $\overline{C}$ and $\underline{C}$\\
				\hline
				1 & No & basic load & $0.06$ & $[19.1;26.3]$\\
				2 & Yes & situation I  & $0.02$& $[19.7;33.8]$  \\
				3 & Yes & situation II  & $0.02$& $[19.4;32.7]$  \\
				4 & No & situation I& $0.18$& $[19.4;33.5]$ \\
				5 & No & situation II& $0.18$& $[19.1;27.5]$ \\
				\hline
			\end{tabular}
		\vspace{-0.25cm}
		\end{table}

	\indent We present the clearing results in Table \ref{clear}. While G1 offers the cheapest upward reserve bid and still owns extra capacity/ramping rate, the SO doesn't clear its entire upward reserve bid, instead the more expensive resources G2 and G3 are cleared. The reason is that the extra upward reserve from G1 can't be delivered in scenarios when the branch outage happens. It can also be observed that there are uniform energy and reserve prices at each bus. In the meantime, the load fluctuation payments can calculated according to equation (\ref{fluctuation payment}): $\Pi^d(1)$=\$23.3, $\Pi^d(2)$=\$91.7, $\Pi^d(3)$=\$-23.5. Note that d3 fluctuation payment is negative because d3's fluctuations will hedge the fluctuations of d1 and d2.

		\begin{table}[H]
		\footnotesize
			\centering
        \vspace{-0.1cm}
        \setlength{\abovecaptionskip}{-0pt}
        \setlength{\belowcaptionskip}{0pt}
            \captionsetup{font=footnotesize,textfont=sc}
			\caption{Clearing Results for the 2-Bus System}
			\label{clear}
			\begin{tabular}{ccccccc}

				\hline
				Generator & $g$ & $r_U$ & $r_D$& $\eta^g$ & $\eta^U$ & $\eta^D$\\
				\hline
				G1 & $8.0$ & $2.4$ & $0.8$& $25.4$ & $2.0$ & $2.0$ \\
				G2 & $17.0$ & $1.0$ & $0.0$& $35.7$ & $5.3$ & $3.7$ \\
				G3 & $0.0$ & $4.0$ & $0.0$& $35.7$ & $5.3$ & $3.7$ \\
				\hline
			\end{tabular}
		\vspace{-0.2cm}
		\end{table}

		\begin{table}
		\footnotesize

            \captionsetup{font=footnotesize,textfont=sc}
            \caption{Money Flow for the 2-Bus System(\$)}
			\label{money flow}
			\centering
			\vspace{-0.2cm}
			\begin{tabular}{ccccccccc}
				\hline
				& Base & S1 & S2 & S3 & S4 & S5& Total\\
				\hline
				$\Gamma^d$ & $446.9$ & $24.9$ & $25.2$ &$9.7$ & $238.2$  &  $86.0$ & $830.7$\\
				$\Pi^d$\ & $0$ &  $0$ & $8.0$ & $0.8$ & $75.9$   & $6.9$& $91.4$\\
				$\epsilon\Phi^d$ & $0$ & $0$ & $2.1$  & $0$ & $6.4$  & $0$& $8.6$\\
				\hline
				$\Gamma^g$ & $443.4$& $20.0$ &  $23.6$  & $9.7$ & $227.5$  & $86.0$& $810.1$\\
				$\Gamma^U$ & $0$ & $0$ & $2.6$ & $0$  & $28.7$  & $0$& $31.3$\\
				$\Gamma^D$& $0$ & $1.6$ & $0$ & $0$ & $0$ & $0$ & $1.6$\\
				$\epsilon\Phi^U$ & $0$ & $1.3$ &$4.0$ &$0.9$ & $38.5$ & $7.5$ & $52.1$\\
				$\epsilon\Phi^D$ & $0$ & $0.9$ & $0.1$ & $0.1$ & $0$ & $0.6$ & $1.7$\\
				$\Delta$ & $3.5$& $2.9$ &$1.0$ &$0$& $12.7$ & $0$ & $20.1$\\
				\hline
			\end{tabular}
			\begin{tablenotes}
			\item[a] base / S1-S5: denote the base case / scenario 1-5;
			\item[b] Revenue adequacy: $(\Gamma^d=\Gamma^g+\Delta)$ holds for column 1, $(\Gamma^d\!+\!\Pi^d\!-\!\epsilon\Phi^d\!=\!\Gamma^g\!+\!\Gamma^U\!+\!\Gamma^D\!+\!\epsilon\Phi^U\!-\!\epsilon\Phi^D\!+\!\Delta)$ holds for other columns.
			\end{tablenotes}
			\vspace{-0.3cm}
		\end{table}

	\indent In Table \ref{money flow}, the money flow is presented. Explicitly revenue adequacy holds in the base case, in each scenario and in total.

	\subsection{IEEE 118-Bus System}
	\indent Simulations on the modified IEEE 118-bus system are also reported. Outage probabilities of lines 21, 55, and 102 are set to be 10\%. The line capacity exceeding rates are set to be 1.2 for all scenarios. The original load 59 will be equally separated into two loads: new load 59 and load 119. There will be two load fluctuation situations, each with occurrence probability 10\%, and the fluctuation levels of all loads will be 3\% in both situations: in situation I, d119 will increase by 3\% while other loads will decrease by 3\%; in situation II, d119 will decrease by 3\% while others will rise by 3\%. The generators' energy and reserve bids are modified to be linear. The total revenue inadequacy in this case is $\$2.26$, which is rather small compared to the expected total cost $\$89476.4$.
	\begin{figure}[!t]

		\centering
		\includegraphics[width=3.9in]{plot_M_on_fluc_d_1_2.pdf}
		\vspace{-0.25cm}
        \setlength{\abovecaptionskip}{-1.5pt}
        \setlength{\belowcaptionskip}{-3.5pt}

	    \captionsetup{font=footnotesize,singlelinecheck=false,}
		\caption{Fluctuation Payments from d59 (blue) and d119 (red) with the Fluctuation Levels of d59 (left) and d119 (right)}
		\label{plot_M_on_fluc_d_1_2}
		\vspace{-0.35cm}
	\end{figure}

     \indent In Fig. \ref{plot_M_on_fluc_d_1_2}, the fluctuation payment from d119 is negative because its fluctuation will hedge the others' fluctuations in all scenarios, so d119's fluctuations will be credited. For the left figure, With the rising fluctuation level of d59, the fluctuation payment from d59 will increase while the fluctuation credit to d119 will increase because d59's rising fluctuation level will bring more uncertainties, therefore enhancing the value of d119's fluctuations hedge against others' fluctuations. For the right figure, while the fluctuation credit to d119 will increase with its rising fluctuation level, fluctuation payment from d59 will decrease because the rising fluctuation level of d119 will reduce the impact of d59's fluctuations on system balance.

		    \begin{figure}[t]
		\centering

		\includegraphics[width=3.9in]{plot_M_on_u_b.pdf}
		\vspace{-0.15cm}

		\vspace{-0.9cm}
		\captionsetup{font=footnotesize,singlelinecheck=false}
		\caption{Price of $r_U(7)$ with the Increasing Fluctuation Levels of d15 (left) and Price of $r_D(28)$ with the Increasing Fluctuation Levels of d66 (right)}
		\label{plot_M_on_u_b}
		\vspace{-0.6cm}
	\end{figure}

    Fig. \ref{plot_M_on_u_b} depicts that with the increasing fluctuation level of d15, the upward reserve price of G7 at bus 15 will increase (left), and with the increasing fluctuation level of d66, the downward reserve price of G28 at bus 66 will increase (right). Underlining in these phenomenons is that the increasing fluctuation levels of d15 and d66 will bring in more uncertainties, therefore enhancing the value of reserve at the same bus.

	\section{Conclusions}
	In this paper, we proposed a scenario-oriented energy-reserve co-optimization model which considers the re-dispatch costs in all non-base scenarios to minimize the expected system total cost, and includes the network constraints in all scenarios to ensure the deliverability of reserve. We defined energy and reserve marginal prices which are locational uniform prices under certain assumptions, and proposed the associated settlement process which can guarantee revenue adequacy of the system operator. In future studies, we aim to include generator outages into the model, and consider the coupling between reserve and ramping in multi-period operations.

	\ifCLASSOPTIONcaptionsoff
	\newpage
	\fi

	%

	%

	%

	\bibliographystyle{ieeetr}


\section{\label{}}

A droplet of a Newtonian liquid which impinges on a highly solvophobic surface, will normally rebound, at least partially.
This rebound is reliably suppressed by the addition of a small amount of polymer, for example, 200 ppm of polyethylene oxide (PEO) in water.\cite{bertola2013}
The anti-rebound behavior of a dilute polymer solution droplet is of great scientific and practical interest, but its underlying physics is poorly understood.\cite{wirth1991, monteux2008}
The fundamental interest arises particularly from the observation that adding very low concentration of polymer to a droplet does not alter its shear viscosity and liquid-vapor surface tension, but leads to a dramatic change in its tendency to rebound.
This contradicts the traditional explanation of spreading and retraction processes of an impacting droplet in terms of a balance between inertia, viscous force, and capillary force.
They are typically characterized by the Reynolds number, Re = $D_{0}v_\mathrm{imp}\rho/\eta_\mathrm{s}$ and the Weber number, We = $D_{0}v_\mathrm{imp}^2\rho/\gamma$, where $D_0$ is the equilibrium droplet diameter, $v_\mathrm{imp}$ is the impact velocity, $\rho$ is the droplet mass density, $\eta_\mathrm{s}$ is the shear viscosity, and $\gamma$ is the fluid's liquid-vapor surface tension.\cite{pasandidehfard1996,josserand2016,ukiwe2005,attane2007,lee2016}
However, there is no concrete relation between these numbers and the rebound tendency of a dilute polymer solution droplet.
This calls for new theoretical and numerical approaches.

Aided by the development of high speed imaging techniques, an early attempt to understand the polymer effect on the rebound tendency was carried out.
It suggested that non-Newtonian elongational viscosity brought by the polymer leads to a reduced retraction velocity of a droplet.\cite{bergeron2000}
This argument, however, has been debated, because the liquid is under shear flow, not under elongational flow, during the retraction stage.
A later study attributed the rebound suppression to the increased friction of the three-phase contact line during the retraction stage, imposed by polymer deposition on the substrate.\cite{zang2013}
This view was also supported by stick-slip dynamics of the contact line and lamella-structured polymer deposit left on the substrate.\cite{smith2010, bertola2013}
An experiment demonstrating that asymmetrical pre-deposition of polymer on a substrate slows the retraction speed, with the consequence of an asymmetric rebound, also highlighted the role of the additional friction provided by the deposited polymer.\cite{huh2015}
Most recently, the importance of the elongational viscosity was revisited by a high damping coefficient of height oscillation of a droplet of PEO in water after impact.\cite{chen2018}
However, the interpretation of the rebound in above studies is limited to the phenomenology and it lacks a molecular picture.
No direct evidence of a polymer effect on the anti-rebound has been provided by numerical solutions of the Navier-Stokes equation with finite-element nonlinear elastic polymer models.\cite{izbassarov2016,tembely2019}
This has been attributed to the inability to describe the wettability change of a substrate and the variation of local polymer concentration within the droplet.
A more suitable way to find the physics behind the anti-rebound effect is a particle-based simulation, but to our best knowledge, no such studies have been reported.
In this letter, we investigate the anti-rebound behavior of an impacting droplet of a dilute polymer solution, while varying systematically the molecular weight and the concentration of the polymer by using multi-body dissipative particle dynamics (MDPD) simulations.

The central idea of the original dissipative particle dynamics (DPD) simulation is to reproduce correct hydrodynamics by viscous forces between groups of atoms.\cite{espanol1995,espanol2017}
In DPD, a large integration time step is possible due to the soft nature of purely repulsive potential functions, thus mesoscopic length and time scales are achievable.
Recently, keeping all the benefits of the DPD method, the limitation of the invariable quadratic equation of state of DPD was overcome by the addition of a local density-dependent multi-body attractive force, which is named MDPD.\cite{warren2003}
The MDPD allows liquid-vapor coexistence and generates a sharp interface, which is necessary for the droplet simulations in this work.
We investigate droplets of six different polymer molecular weights (number of polymer beads $N_\mathrm{p}$= 1--100), and of five different number fractions of polymer beads ($x_\mathrm{p}$= 0.5--2.5 $\times 10^{-2}$) for each $N_\mathrm{p}$.
The highest polymer concentration is only one tenth of the overlap concentration even for $N_\mathrm{p}$= 100.
Thus, our solutions are in the very dilute regime.
Each droplet contains 4$\times 10^4$ beads including polymer and solvent beads.
The substrate (surface) is modeled by frozen particles on a square lattice.
The repulsion amplitude and the cutoff distance are chosen to be same for all particle pairs ($B=25$, $r_B=1.0$). (See Supplemental Material for details of simulation methods and models, including DPD reduced units)
The amplitudes of the attractive force ($A$) differ depending on the pair type such that $A_\mathrm{s/s}=A_\mathrm{p/p}=A_\mathrm{s/p}=-40$, $A_\mathrm{s/w}=-10$, and $A_\mathrm{p/w}=-30$, where subscripts s, p, w refer to solvent, polymer and substrate (``wall'') beads, respectively.
Such a choice of interaction parameters at temperature $k_\mathrm{B}T=1$ mimics polymer in a good solvent, and a solvophobic but polymerphilic substrate.
The equilibrium contact angle for a pure solvent droplet on the substrate is 155\textdegree.
After equilibrating a free droplet, we let it impact on the substrate with a constant velocity of 2.
The surface tension of free droplets and the zero-shear viscosity of corresponding polymer solutions increase slightly as either $N_\mathrm{p}$ or $x_\mathrm{p}$ increases, but without meaningful difference beyond the statistical error.
Together with equilibrium droplet diameters ($D_{0}=23.2$) and liquid mass densities ($\rho=6.10$) of droplets, which are nearly constant irrespective of polymer composition, the Reynolds and Weber numbers are calculated to be similar as Re= 35--42 and We= 78--81. (See Supplemental Material)

\begin{figure*}[t]
\centering
\includegraphics[width=1.0\textwidth]{outcomes.eps}
\caption{\label{fig-outcomes}Simulation snapshots of four different outcomes from side view: (a) complete rebound, (b) partial rebound, (c) deposition, and (d) deposition by polymer necking. In all images, solvent and substrate particles are colored in cyan and yellow, respectively. Each polymer strand is depicted in different color. In Figure (e), the most probable outcomes are shown for different polymer molecular weight ($N_\mathrm{p}$) and concentration ($x_\mathrm{p}$). Two outcomes for a given $N_\mathrm{p}$ and $x_\mathrm{p}$ (for example, deposition and partial rebound for $N_\mathrm{p}$= 20, $x_\mathrm{p}$= 0.015) indicate that two outcomes are equally probable out of five independent trajectories. The black dashed line on the graph is a rough boundary between rebound and deposition. Snapshots are obtained by VMD.\cite{humphrey1996}}
\end{figure*}

Right after impact the droplet first spreads.
After reaching a maximum expansion, it retracts back toward the droplet center followed by four different outcomes of rebound: complete rebound, partial rebound, deposition, and deposition by polymer necking as displayed in Figures \ref{fig-outcomes}(a)--(d). (See also movies in Supplemental Material)
Figure \ref{fig-outcomes}(e) shows the most probable outcome out of five independent trajectories for each polymer composition.
A droplet of pure solvent rebounds completely resulting from high enough inertia and surface tension as intended by chosen simulation parameters.
Droplets of small $x_\mathrm{p}$ and small $N_\mathrm{p}$ (left lower part of Figure \ref{fig-outcomes}(e)) partially rebound with a small satellite droplet left on the substrate, which was also observed in several experiments.\cite{bertola2013,huh2015}
The satellite drop is enriched in polymer of the number fraction larger than 0.5.
As in experiments, we find a general trend that droplets containing polymer of longer $N_\mathrm{p}$ and higher $x_\mathrm{p}$ do not rebound, although their solution viscosity is not much larger than that of the pure solvent.
A key factor of rebound suppression, which can be extracted from this figure is that polymer molecules are largely adsorbed on the substrate to mediate the unfavorable interaction between substrate and solvent.
Moreover, the remaining part of the polymer not adsorbed on the substrate is stretched normal to the substrate, i.e. the direction of rebound, with the flowing solvent.
After the polymer is highly stretched, it tends to recover its coiled structure, exerting a force against the rebound to the surrounding solvent particles.
The effect of the stretched polymer is most prominent in the final outcome, namely deposition by polymer necking (Figure 1(d)).
This is observed especially for the longest polymer of $N_\mathrm{p}=$ 50 or 100 at an intermediate concentration of $x_\mathrm{p}=$ 0.010 or 0.015.
As only some polymer beads are adsorbed on the substrate, the majority of a long polymer strand can stretch in the direction of rebound.
Most solvent particles temporarily leave the substrate and form a partial droplet above the substrate.
The polymer, however, whose one end is adsorbed on the substrate and the other end is buried in the solvent droplet, actually pulls the droplet back while recovering its coiled structure ($t=$ 200).
As a result, the whole droplet ends up deposited on the substrate.
Such a polymer necking conformation has not been found or postulated in any droplet impact experiments, but a droplet shape very similar to our simulation snapshot at $t=$ 200 in Figure \ref{fig-outcomes}(d) has been observed in experiments using a high-speed camera.\cite{huh2015}
A stretched polymer configuration and the necking morphology observed in droplet pinch-off experiments\cite{sachdev2016} also strongly support the reality of the formation of a polymer necking morphology during droplet impact.

\begin{figure}[t!]
\centering
\subfigure[]{\label{fig-spreading}\includegraphics[width=0.48\textwidth]{spreading_factor.eps}}
\subfigure[]{\label{fig-velocity}\includegraphics[width=0.48\textwidth]{velocity.eps}}

\caption{\label{fig-sprvel} (a) Spreading factor and  (b) average particle axial velocity as a function of time for droplets of different outcomes: complete rebound (R), partial rebound (P), deposition (D) and deposition by necking (N). Vertical dashed lines indicate the beginning of a hopping stage at $t=60$. Arrows in Figure (b) indicate the hopping velocity.}
\end{figure}

The spreading diameter $D$ is the diameter of the circle formed by the three-phase contact line.
It is shown as the spreading factor, $\beta=D/D_0$ as a function of time in Figure \ref{fig-spreading}.
The maximum spreading factors, $\beta_\mathrm{max}$, of all droplets are almost the same independent of polymer composition because they have almost the same We and Re numbers.\citep{ukiwe2005}
This is also evident in the snapshots at $t=$ 27 of Figure \ref{fig-outcomes}(a)--(d).
Figure \ref{fig-spreading} shows that the velocity of the three-phase contact line is reduced for non-rebounding droplets after $\beta_\mathrm{max}$ is attained.
However, there is a striking difference between our simulation and experimental results: In the simulation, the slowing of the contact line motion begins at the very end of the retraction stage, $t\approx 60$ (vertical dashed lines in Figure \ref{fig-sprvel}) when $\beta\approx 1$, while in experiments the retraction slows down right after $\beta$ reaches $\beta_\mathrm{max}$.\cite{bergeron2000, bertola2013, huh2015}
This result is consistent with the fact that droplets, which later will have different rebounding outcomes, still show almost the same morphology and size at early times up to $t=$ 60 in Figures \ref{fig-outcomes}(a)--(d).
This behavior is more clearly seen in the average particle velocity in the direction normal to the surface in Figure \ref{fig-velocity}.
This figure shows that particles in the non-rebounding droplet also reach an axial velocity as high as in the rebounding droplets at $t\approx 60$.
Only after that time, the velocity decreases depending on the strength of resistance against the rebound motion.
Looking at the shape of the droplets at $t=$ 60 in Figure \ref{fig-outcomes}(a)--(d), they are no longer compressible in the surface plane and the particles' momentum starts to be diverted from parallel to perpendicular direction to the substrate. (See also Supplemental Material for almost zero radial particle velocities at $t>60$)
Therefore, the rebound of the droplets is suppressed not by the resistance against the parallel retraction but by the resistance against the perpendicular motion.
We newly name this stage of $t>60$ as a \textit{hopping stage} and will focus on this stage in the rest of this letter.
One should keep in mind that our finding does not deny the experimental observations but rather states that the rebound can also be suppressed during the hopping stage, even though the retraction velocity is the same as for the rebounding droplet.

\begin{figure}[t!]
\subfigure[]{\label{fig-normalstress}\includegraphics[width=0.48\textwidth]{normalstress.eps}}
\subfigure[]{\label{fig-contangle}\includegraphics[width=0.48\textwidth]{contactangle.eps}}
\subfigure[]{\label{fig-polads}\includegraphics[width=0.48\textwidth]{poladsorption.eps}}
\caption{\label{fig-polcont} (a) First normal stress difference, (b) contact angle, and (c) the number of adsorbed polymer beads as a function of time for droplets of different outcomes. Vertical dashed lines indicate the beginning of a hopping regime at $t=60$. 5-point running averages of the contact angle for better statistics are plotted in (b). Horizontal black dotted lines in (a) and (b) indicate 0 and the equilibrium contact angle of the pure-solvent droplet, $\theta_\mathrm{eq}=155$, respectively.}
\end{figure}

To monitor the force that droplet beads feel along the perpendicular direction to the substrate ($z$-direction), we calculate the first normal stress difference, $N_1=\sigma_{zz}-(\sigma_{xx}+\sigma_{yy})/2$ in Figure \ref{fig-normalstress}.
This Figure shows $N_1$ obtained from 5-point running averages of stress tensor components.
Consistent with the finding in Figure \ref{fig-sprvel}, the main difference in $N_1$ between different rebounding outcomes appears during the hopping stage.
The negative stress in this stage represents a net attractive force between particles and its amount is bigger for non-rebounding droplets than for rebounding droplets especially at large $t>60$.
As the elongational viscosity, $\eta_\mathrm{e}$, of a liquid under uniaxial elongational flow is given by $\eta_\mathrm{e}=-N_1 / \dot{\epsilon}$, where $\dot{\epsilon}$ is an elongational flow rate, the negative $N_1$ is mainly contributed by the elongational force of polymer.
However, even the pure-solvent droplet has a negative stress at the early hopping stage because the net attractive force primarily stems from interaction at the liquid-solid interface.(Solvent and substrate beads also weakly attract each other.)
At the interface, if polymer beads are adsorbed on the substrate as shown in Figure \ref{fig-outcomes}, they attract both solvent and substrate beads to mediate unfavorable solvent-substrate interaction, leading to large $|N_1|$.

A time evolution of the contact angle shows the effect of deposited polymer on the dynamic wetting. (See Supplemental Material for the method)
We find $\approx$10\textdegree~ difference between rebounding and non-rebounding droplets in the retracting (receding) contact angle at $t=60$, which is similar to the experimental observation.\cite{bertola2015}
However, a big difference is again shown in the hopping stage.
While the contact angle of the completely rebounding droplet stays $\approx$ 155\textdegree~ until it detaches, the contact angle of non-rebounding droplets decreases to $\approx$ 90\textdegree~ until $t\approx 100 $ which is an equilibrium contact angle of a droplet containing only polymer beads.
The decreased contact angle during the hopping stage means strong attractive forces at the liquid-solid interface.
Taking an advantage of particle-based simulations, we calculate the number of polymer beads adsorbed on the substrate ($n_\mathrm{ads}$) as a function of time (Figure \ref{fig-polads}), which quantifies the change of the wettability of the substrate.
We define as adsorbed polymer beads those, whose distance from the substrate is less than 1, the cutoff distance of the attractive force.
Previous numerical studies showed that there is no universal behavior of polymer adsorption and desorption under shear flow, and both depend strongly on the flow rate and the polymer-substrate interaction.\cite{dutta2013, dutta2015}
Under our specific condition $n_\mathrm{ads}$ increases during the spreading and early retraction stages, and stays constant at the late retraction stage ($t=50-60$) due to the balance between flow rate and the substrate-polymer attraction.
At the beginning of the hopping stage ($t=60)$) the liquid flow direction changes to the substrate normal, and $n_\mathrm{ads}$ starts to decrease.
Therefore, we take the maximum value of $n_\mathrm{ads}$, $n_\mathrm{ads,max}$, to characterize the change of the substrate wettability by the polymer adsorption.
As expected, we obtain larger $n_\mathrm{ads,max}$ for longer $N_\mathrm{p}$ and larger $x_\mathrm{p}$, (See Supplemental Material) which represents the stronger attractive force acting on a droplet during the hopping stage.

\begin{figure}[t]
\centering
\includegraphics[width=0.48\textwidth]{hopvel_vs_pads.eps}
\caption{\label{fig-hopvel} Hopping velocity versus the maximum value of the number of adsorbed polymer beads for different $N_\mathrm{p}$ and $x_\mathrm{p}$. Each $v_\mathrm{hop}$ is obtained from each trajectory in which the droplet rebounds. Fitted linear functions for each $N_\mathrm{p}$ ($v_\mathrm{hop}=-sn_\mathrm{ads,max}+c$) are represented by dashed lines, where the fitted $s$ are given in the legend. }
\end{figure}

To quantitatively understand the effect of polymer adsorption on the rebound behavior, we plot the hopping velocity, $v_\mathrm{hop}$, as a function of $n_\mathrm{ads,max}$ in Figure \ref{fig-hopvel}.
Here, $v_\mathrm{hop}$ is defined as the constant velocity after the droplet detaches from the substrate (Figure \ref{fig-velocity}).
As expected, $v_\mathrm{hop}$ decreases as more polymer beads are adsorbed on the substrate for a given $N_\mathrm{p}$.
We can fit the data for each $N_\mathrm{p}$ with a line (dashed lines), especially for short polymers.
It is noticeable that the slope monotonously decreases with increasing $N_\mathrm{p}$.
In other words, for a given $n_\mathrm{ads,max}$ characterizing a given substrate wettability, a droplet containing longer polymer chains is less likely to rebound.
This finding clearly indicates that the polymer elongation force additionally affects the rebound behavior, as already expected from Figures \ref{fig-outcomes} and \ref{fig-normalstress}.
The elongation force originates from the tendency of the polymer to recover its coiled equilibrium conformation, and it generates a pulling force on the droplet by hydrodynamic interaction.
Assuming Hookean elasticity of the stretched polymer in infinite dilution, the recovering force is proportional to the end-to-end distance or alternatively the stretching factor of the stretched polymer.
We find that the stretching factor $\langle S_z \rangle$, which is the difference along the $z$-direction between uppermost and lowermost beads of a chain, is roughly proportional to $N_\mathrm{p}$ because the polymer is almost fully stretched. (See Supplemental Material)
This indicates that the elongation force per chain also scales as $N_\mathrm{p}$.
However, the number of polymers for a given $x_\mathrm{p}$ is inversely proportional to $N_\mathrm{p}$, leading to the same total elongation force under assumptions of no interchain interactions and a homogeneous polymer concentration.
The main difference in the resistance comes from how long the elongation force acts on the droplet which is again related to the amount of stretching.
This can be also observed in Figure \ref{fig-normalstress} where the minimum value of the negative stress during the hopping stage is almost the same for all non-rebounding droplets and the main difference is how long the negative stress lasts.
As a result, a higher polymer molecular weight leads to a stronger total resistance against rebound.
The slope of the fitted functions ($s$) in Figure \ref{fig-hopvel} scales as $\sim -N_\mathrm{p}^{0.17(3)}$. (See Supplemental Material for a graph of $v_\mathrm{hop}$ as a function of $n_\mathrm{ads,max}\cdot N_\mathrm{p}^{0.17}$)
However, a droplet containing longer polymers has a larger polymer concentration locally near the neck than of shorter polymer,  which makes it even more complicated to find a relation.
We are now investigating the rebound behavior depending on the elongational viscosity of the polymer solutions by considering the local polymer concentration and the stretching of the unadsorbed polymer strand during the hopping stage.

To conclude, we studied the rebound of droplets of polymer solution from a solvophobic substrate by MDPD simulations.
We found that the droplet rebound is suppressed by polymer after the end of the retraction stage when the main flux has been converted from parallel to perpendicular to the substrate.
While the velocity of the three-phase contact line during the retraction stage is not reduced, the polymer manifests itself during this hopping stage.
It makes the two polymer contributions to reducing the rebound tendency: The first is to alter the liquid-solid interfacial interaction by adsorbed polymer primarily generating an attraction.
The second is the restoring force of the polymer after it has been stretched in perpendicular direction.
The alteration of the liquid-solid interface is characterized by the number of polymer beads adsorbed on the substrate that shows an anti-correlation with the hopping velocity.
For a given substrate wettability modified by polymer adsorption, a long polymer pulls the droplet back to the substrate for a long time because it is highly stretched.
Thus, we find that the polymer counteracts droplet rebound at the stage, when the droplet wants to detach from the substrate, rather than early, when the droplet retracts from its maximum lateral expansion.
While our interpretation is consistent with the fact that a small polymer content does not significantly changes the viscosity and the surface tension of the liquid, we cannot rule out that other mechanisms might prevail for different combinations of material's properties, such as solvophilicities vs. adsorptivities.
At the same time, our results open the way to tailoring polymer additives against droplet rebound.

\begin{acknowledgements}

We would like to thank for the financial support by the Collaborative Research Centre Transregio 75 (project number 84292822)
\textit{Droplet Dynamics in Extreme Environments} of the Deutsche Forschungsgemeinschaft. Calculations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt and the GCS Supercomputer JUWELS at J\"ulich Supercomputing Centre (JSC) through the John von Neumann Institute for Computing (NIC).

\end{acknowledgements}