text:
\begin{document}
\begin{abstract}
We characterize Zagier's generating series of extended period polynomials of normalized Hecke eigenforms for $\PSL_2(\bZ)$ in terms of the period relations and existence of a suitable factorization. For this we use a characterization of the Kronecker function as the ``fundamental solution'' of the Fay identity.
\end{abstract}
\maketitle

\section{Introduction}

\subsection{Zagier's generating series of extended period polynomials}

To every cusp form\footnote{Here and in what follows every modular form will be for the full modular group $\PSL_2(\bZ):=\SL_2(\bZ)/\{\pm 1\}$.} $f$ one can attach its period polynomial $r_f(X)$. This is an important arithmetic invariant whose coefficients are, up to elementary factors, the critical values of the completed L-function of $f$. Period polynomials were studied extensively by Eichler, Shimura and Manin (see for example \cite[Chapter V]{Lang:ModularForms} for a textbook account) and their definition was extended to all modular forms by Zagier, \cite{Zagier:Periods}. Nowadays, extended period polynomials are well-studied objects in number theory which keep attracting attention. For example, they are related to zeta functions of real quadratic fields, \cite{KZ}, and multiple zeta values, \cite{GKZ}. More recently, an analogue of the Riemann hypothesis was proved for them, \cite{CFI,Ono}.

In \cite{Zagier:Periods}, Zagier introduced a generating series $C(X,Y;\tau;T)$ of extended period polynomials of normalized Hecke eigenforms (see Section \ref{ssec:4.1} for its definition), and related it to the Kronecker function, \cite[Chapter VIII]{Weil:Elliptic},
\begin{equation}
F_{\tau}(u,v):=\frac{\theta'_{\tau}(0)\theta_{\tau}(u+v)}{\theta_{\tau}(u)\theta_{\tau}(v)}.
\end{equation}
Here $\theta_{\tau}(u)$ denotes Jacobi's odd theta function (see Section \ref{ssec:2.1} for the definition). More precisely, Zagier's main result is as follows.
\begin{thm}[Zagier; {\cite{Zagier:Periods}, eqn. (17)}] \label{thm:Zagier}
We have the following identity of formal series in variables $X,Y,T$:
\begin{equation} \label{eqn:ZagiersTheorem}
C(X,Y;\tau;T)=F_{\tau}(T,-XYT)F_{\tau}(XT,YT).
\end{equation}
\end{thm}
As two applications, we mention the algorithmic computation (up to scalars) of the period polynomials of all cuspidal Hecke eigenforms, as well as a new proof of the Eichler--Selberg trace formula for $\PSL_2(\bZ)$, \cite{Zagier:Hecke}.
\subsection{Period relations and the main result}
The goal of this paper is to give another application of Zagier's result, namely a purely algebraic characterization of $C$ in terms of the period relations for extended period polynomials
\begin{equation} \label{eqn:periodrels}
\begin{aligned}
r_f(X)\vert(1+S)=r_f(X)\vert(1+U+U^2)=0,
\end{aligned}
\end{equation}
where $S=\left(\begin{smallmatrix}0&-1\\1&0\end{smallmatrix}\right)$, $U=\left(\begin{smallmatrix}1&-1\\1&0\end{smallmatrix}\right)$, and $\vert$ denotes the slash operator. Indeed, since $C$ is defined as a generating series of extended period polynomials it also satisfies a suitable version of the period relations, see Equation \eqref{eqn:periodrelations}.

It turns out that the period relations alone do not suffice to characterize $C$, as for any modular form $f$ of weight $k$ the series $C+r_f(X)T^{k-2}$ also satisfies the period relations. However, if we demand in addition the existence of a factorization as in Theorem \ref{thm:Zagier}, then the situation is much better and our main result is as follows.
\begin{thm} \label{thm:main2}
For a formal Laurent series $f(u,v)=\sum_{m,n>\!\!>-\infty}a_{m,n}u^mv^n \in \bC((u,v))$, assume that
\begin{equation}
C_f(X,Y,T):=f(T,-XYT)f(XT,YT) \in \bC((X,Y,T))
\end{equation}
satisfies the period relations, Equation \eqref{eqn:periodrelations}. Then either:
\begin{enumerate}
\item[(i)]
There exist unique $\alpha,\beta \in \bC$ such that
\begin{equation} \label{eqn:case1}
C_f(X,Y,T)=C_{P_{\alpha,\beta}}(X,Y,T)=\frac{(\alpha XY-\beta)(\beta X+\alpha Y)}{X^2Y^2T^2},
\end{equation}
where $P_{\alpha,\beta}(u,v)=\alpha u^{-1}+\beta v^{-1}$.
\item[(ii)]
There exist $\alpha,\beta,\delta \in \bC^{\times}$ and $\tau \in \fH \cup \{i\infty\}$ such that
\begin{equation} \label{eqn:case2}
C_f(X,Y,T)=C_{F_{\tau}}(X,Y',T')=C(X,Y';\tau;T'), \quad \mbox{where } T':=\frac{T}{\alpha}, \, Y':=\frac{\alpha Y}{\beta},
\end{equation}
Moreover, if $\alpha',\beta',\delta',\tau'$ are a different choice of parameters as above, then there exists a matrix $\left(\begin{smallmatrix}a&b\\c&d\end{smallmatrix}\right) \in \SL_2(\bZ)$ such that
\begin{equation}
(\alpha',\beta',\delta',\tau')=\left( \alpha(c\tau+d),\beta(c\tau+d),\frac{\delta}{c\tau+d},\frac{a\tau+b}{c\tau+d} \right).
\end{equation}
\end{enumerate}
\end{thm}
The basic idea is that, via Theorem \ref{thm:Zagier} above, the period relations for $C$ are equivalent to the Fay identity for $F_{\tau}(u,v)$, Equation \eqref{eqn:Fay}, which is a consequence of Fay's trisecant identity for theta functions, \cite[p.34, eqn. (45)]{Fay:Theta}. More precisely,  $C_f$ satisfies the period relations if and only if $f$ satisfies the Fay identity, see Proposition \ref{prop:periodsFay}.
We are therefore led to characterize solutions of the Fay identity, for which we have the following result.
\begin{thm} \label{thm:main}
Assume that $f(u,v)=\sum_{m,n>\!\!>-\infty}a_{m,n}u^mv^n$ satisfies the Fay identity, Equation \eqref{eqn:Fay}.
\begin{enumerate}
\item[(i)]
If $(a_{3,0},a_{5,0})=(0,0)$, then there exist unique $\alpha,\beta,\gamma \in \bC$ such that
\begin{equation} \label{eqn:thmmainpolar}
f(u,v)=e^{\gamma uv}\left(\frac{\alpha}{u}+\frac{\beta}{v}\right).
\end{equation}
\item[(ii)]
If $(a_{3,0},a_{5,0}) \neq (0,0)$, then there exist $\alpha,\beta,\delta \in \bC^{\times}$, $\gamma \in \bC$ and $\tau \in \fH \cup \{i\infty\}$ such that
\begin{equation} \label{eqn:thmmainelliptic}
f(u,v)=\delta e^{\gamma uv}F_{\tau}\left(\frac u\alpha, \frac v\beta\right).
\end{equation}
Moreover, if $\alpha',\beta',\gamma',\delta',\tau'$ are a different choice of parameters as above, then there exists $\left(\begin{smallmatrix}a&b\\c&d\end{smallmatrix}\right) \in \SL_2(\bZ)$ such that
\begin{equation}
(\alpha',\beta',\gamma',\delta',\tau')=\left( \alpha(c\tau+d),\beta(c\tau+d),\gamma-\frac{c/2\pi i}{\alpha\beta(c\tau+d)},\frac{\delta}{c\tau+d},\frac{a\tau+b}{c\tau+d} \right).
\end{equation}
\end{enumerate}
\end{thm}
The idea of proof is that the Fay identity implies recurrence relations for the coefficients of $f$ such that every solution is uniquely determined by at most five of its coefficients (see Proposition \ref{prop:Faybasic2} and Theorem \ref{thm:generators} for the precise statements). Varying the parameters $\alpha,\beta,\gamma,\delta,\tau$ suitably, we can arrange that $f$ is equal to one of the functions in Theorem \ref{thm:main}.

After a previous version of this manuscript had been submitted, the author was informed that the existence part of Theorem \ref{thm:main} is equivalent to a result of Polishchuk, \cite[Theorem 5]{Polishchuk:YB}. To be precise, Polishchuk works with the (scalar) associative Yang--Baxter equation which is equivalent to the Fay identity, \cite{Polishchuk:Tangents}, and classifies solutions in the space of meromorphic functions defined in a neighborhood of $(0,0)$ instead of formal Laurent series.
The two proofs are quite similar but ours may still be of independent interest as it yields slightly more information about the coefficients of solutions to the Fay identity (see Proposition \ref{prop:coefficients}). This in turn sheds some light on how solutions to the associative Yang--Baxter equation can be constructed algorithmically.

\subsection{Content}
Section \ref{sec:2} is preliminary; we recall the definition of the Kronecker function $F_{\tau}(u,v)$ and state the Fay identity. This identity is then studied in detail in Section \ref{sec:3} which culminates in the proof of Theorem \ref{thm:main}. In Section 4 we explain the relation between Fay identity and period relations, building on Zagier's fundamental result, and finish by giving a proof of Theorem \ref{thm:main2}.
\smallskip

{\bf Acknowledgements:} This work was done while the author was a JSPS Postdoctoral Fellow, partly supported by JSPS KAKENHI Grant No. 17F17020. The author would like to thank his academic host, Professor Masanobu Kaneko, for his support and fruitful discussions on the contents of this manuscript, as well as Ulf K\"uhn and the referees for crucial comments and helpful feedback. Final corrections were made while the author was a Postdoctoral Research Assistant at University of Oxford, partly supported by ERC grant 724638.
\smallskip

{\bf Notation:}
Given variables $X_1,\ldots,X_n$ and a ring $R$, we will denote by $R((X_1,\ldots,X_n))$ the $R$-algebra of formal Laurent series of the form
\begin{equation}
\sum_{m_1,\ldots,m_n>\!\!>-\infty}a_{m_1,\ldots,m_n}X_1^{m_1}\ldots X_n^{m_n}, \quad \mbox{with } a_{m_1,\ldots,m_n} \in R.
\end{equation}

\section{The Kronecker function and its functional equations} \label{sec:2}

The reference for this section is \cite{Zagier:Periods}.

\subsection{The Fay identity for the Kronecker function} \label{ssec:2.1}
Consider the classical odd Jacobi theta function $\theta_{\tau}(u):=\sum_{n\in \bZ}(-1)^nq^{\frac 12\left(n+\frac 12\right)^2}e^{\left(n+\frac 12\right)u}$. It is entire, has simple zeros exactly for $u \in 2\pi i(\bZ+\bZ\tau)$ and satisfies $\theta_{\tau}(u+2\pi i(m\tau+r))=(-1)^{m+r}q^{-\frac{m^2}{2}}e^{-mu}\theta_{\tau}(u)$, for $m,r\in \bZ$.
\begin{dfn}
The Kronecker function is the meromorphic function $F_{\tau}: \bC^2 \rightarrow \bC$, defined by
\begin{equation}
F_{\tau}(u,v)=\frac{\theta'_{\tau}(0)\theta_{\tau}(u+v)}{\theta_{\tau}(u)\theta_{\tau}(v)}.
\end{equation}
\end{dfn}
From the properties of $\theta_{\tau}$ mentioned above, it follows that $F_{\tau}(u,v)$ has simple poles if $u$ or $v$ is in $2\pi i(\bZ+\bZ\tau)$, simple zeros for $u+v \in 2\pi i(\bZ+\bZ\tau)$ and satisfies $F_{\tau}(u+2\pi i(n\tau+s),v+2\pi i(m\tau+r))=q^{-mn}e^{-mu-nv}F_{\tau}(u,v)$, for $m,n,r,s\in \bZ$. Also, at the cusp $i\infty$, the Kronecker function degenerates to a trigonometric function:
\begin{equation}
F_{\tau}(u,v)\vert_{\tau=i\infty}=\frac 12\left(\coth\left( \frac u2\right)+\coth\left( \frac v2 \right)\right).
\end{equation}
\begin{rmk}
The above version of the Kronecker function is the one given in \cite{Zagier:Periods}. Other sources such as \cite{BL:MEP,Polishchuk:YB} use the function $2\pi iF_{\tau}(2\pi iu,2\pi iv)$ instead.
\end{rmk}
We are interested in functional equations satisfied by $F_{\tau}$. Since $\theta_{\tau}(-u)=-\theta_{\tau}(u)$, the Kronecker function is antisymmetric, i.e.
$
F_{\tau}(u,v)+F_{\tau}(-u,-v)=0.
$
More interestingly, it satisfies the following three-term functional equation, \cite{BL:MEP}, Proposition 5.(iii).\footnote{Note that our version is slightly different, see Remark \ref{rmk:Fay} for more details.}
\begin{prop} \label{prop:Fay}
The Kronecker function satisfies the Fay identity,
\begin{equation} \label{eqn:Fay}
F_{\tau}(u_1,v_1)F_{\tau}(u_2,v_2)+F_{\tau}(-u_2,v_1-v_2)F_{\tau}(u_1+u_2,v_1)+F_{\tau}(-u_1-u_2,-v_2)F_{\tau}(u_1,v_1-v_2)=0.
\end{equation}
\end{prop}
\begin{proof}
We give a proof for convenience of the reader. Writing out the definition of $F_{\tau}(u,v)$ and multiplying by the common denominator of the left hand side of \eqref{eqn:Fay}, we see that \eqref{eqn:Fay} is equivalent to
\begin{equation}\label{eqn:propFay2}
\begin{aligned}
0&=\theta_{\tau}(u_1+v_1)\theta_{\tau}(u_2+v_2)\theta_{\tau}(u_1+u_2)\theta_{\tau}(v_1-v_2)\\
&-\theta_{\tau}(-u_2+v_1-v_2)\theta_{\tau}(u_1+u_2+v_1)\theta_{\tau}(v_2)\theta_{\tau}(u_1)\\
&+\theta_{\tau}(-u_1-u_2-v_2)\theta_{\tau}(u_1+v_1-v_2)\theta_{\tau}(u_2)\theta_{\tau}(v_1),
\end{aligned}
\end{equation}
where we also used $\theta_{\tau}(-u)=-\theta_{\tau}(u)$. Now substituting
\begin{alignat}{3}
&\alpha_0=u_1+u_2+v_2, \quad &&\alpha_1=-u_2+v_1-v_2, \quad &&\alpha_2=-\alpha_0-\alpha_1=-u_1-v_1,\\
&\beta_0=u_2, \quad &&\beta_1=v_2, &&\beta_2=-\beta_0-\beta_1=-u_2-v_2,
\end{alignat}
and again using antisymmetry of $\theta_{\tau}$, we can write \eqref{eqn:propFay2} in the more symmetric form
\begin{equation}
\sum_{i \in \bZ/3\bZ}\theta_{\tau}(\alpha_i)\theta_{\tau}(\beta_i)\theta_{\tau}(\alpha_{i-1}+\beta_{i+1})\theta_{\tau}(\alpha_{i+1}-\beta_{i-1})=0,
\end{equation}
and this is precisely \cite[Proposition 5]{Zagier:Periods}.
\end{proof}
Slightly more generally, we have the following result.
\begin{cor} \label{cor:solution}
The functions
\begin{equation}
e^{\gamma uv}\left( \frac{\alpha}{u}+\frac{\beta}{v} \right), \quad \mbox{and} \quad \delta e^{\gamma uv}F_{\tau}\left( \frac{u}{\alpha},\frac{v}{\beta} \right),
\end{equation}
both satisfy the Fay identity, for all $\alpha,\beta,\gamma \in \bC$ in the first case and for all  $\alpha,\beta \in \bC^{\times}$, $\gamma,\delta \in \bC$ and $\tau \in \fH \cup \{i\infty\}$ in the second.
\end{cor}
\begin{proof}
The key observation is that the product $e^{\gamma u_1v_1}e^{\gamma u_2v_2}$ is invariant under the following two linear transformations which occur in the Fay identity:
\begin{equation}
\begin{aligned}
(u_1,v_1,u_2,v_2) \mapsto (-u_2,v_1-v_2,u_1+u_2,v_1),\\
(u_1,v_1,u_2,v_2) \mapsto (-u_1-u_2,-v_2,u_1,v_1-v_2).
\end{aligned}
\end{equation}
The corollary then follows from partial fractions in the first case and from Proposition \ref{prop:Fay} in the second.
\end{proof}

\section{Algebraic structure of the Fay identity} \label{sec:3}

In this section, we always let $f(u,v)=\sum_{m,n>\!\!>-\infty}a_{m,n}u^mv^n \in \bC((u,v))$ be a formal Laurent series. The goal of this section is to derive constraints on the coefficients $a_{m,n}$ imposed by the Fay identity. Our main results (Theorem \ref{thm:generators} and Proposition \ref{prop:Faybasic2} below) show that if $f$ satisfies the Fay identity, it is uniquely determined by at most five of its coefficients.
\subsection{Some basic implications of the Fay identity}
We begin by showing how the Fay identity implies the vanishing of many of the coefficients $a_{m,n}$.
\begin{prop} \label{prop:Faybasic}
If $f$ satisfies the Fay identity, then
\begin{equation}
f(u,v)=\frac{a_{-1,0}}{u}+\frac{a_{0,-1}}{v}+\sum_{m,n\geq 0}a_{m,n}u^mv^n, \quad \mbox{with }a_{m,n}=0, \, \mbox{if } m+n \in 2\bZ.
\end{equation}
In particular, $f(-u,-v)=-f(u,v)$. Moreover, if $(a_{-1,0},a_{0,-1})=(0,0)$, then $f=0$.
\end{prop}
\begin{proof}
We first show that if $f \neq 0$, then $f$ must have a pole at either $u=0$ or $v=0$. Indeed, if $f$ does not have a pole at $v=0$, then the Fay identity implies
\begin{align}
f(u_1,0)f(u_2,0)+f(-u_2,0)f(u_1+u_2,0)+f(-u_1-u_2,0)f(u_1,0)&=0,\\
f(u_1,v_1)f(u_2,v_1)+f(-u_2,0)f(u_1+u_2,v_1)+f(-u_1-u_2,-v_1)f(u_1,0)&=0.
\end{align}
The first equation implies $f(u,0)=a_{-1,0}u^{-1}$, as can be seen by comparing coefficients. In particular, if $f$ does not have a pole at $u=0$, we must have $f(u,0)=0$, but then the second equation implies $f(u,v)=0$.

Now assume that $f \neq 0$ and let $M,N \in \bZ$ be the largest integers such that $a_{-M,n} \neq 0$ for some $n \in \bZ$ and $a_{m,-N} \neq 0$ for some $m \in \bZ$. Since $f$ has a pole at either $u=0$ or $v=0$, we have $M\geq 1$ or $N\geq 1$. Define
\begin{equation}
g(v):=(u^Mf(u,v))\vert_{u=0}=\sum_{n>\!\!>-\infty}a_{-M,n}v^n, \quad h(u):=(v^Nf(u,v))\vert_{v=0}=\sum_{m>\!\!>-\infty}a_{m,-N}u^m.
\end{equation}
Both $g$ and $h$ are well-defined by construction of $M$ and $N$. Now, if $f$ has a pole at $u=0$, we multiply the Fay identity by $(u_1u_2)^M$ and then set $u=u_1=u_2=0$ to get
\begin{equation}
g(v_1)g(v_2)+\left(\frac{-1}{2}\right)^Mg(v_1-v_2)g(v_1)+\left(\frac{-1}{2}\right)^Mg(-v_2)g(v_1-v_2)=0,
\end{equation}
and it is straightforward to verify that this implies $M=1$ and that $g(v)=a_{-1,0}$. Likewise, if $f(u,v)$ has a pole at $v=0$ (i.e. $N\geq 1$), then a similar argument yields $N=1$ and $h(u)=a_{0,-1}$. This shows that
\begin{equation}
f(u,v)=\frac{a_{-1,0}}{u}+\frac{a_{0,-1}}{v}+\sum_{m,n \geq 0}a_{m,n}u^mv^n,
\end{equation}
with $(a_{-1,0},a_{0,-1}) \neq (0,0)$ if $f\neq 0$.

It remains to prove that $a_{m,n}=0$ if $m+n$ is even which is equivalent to antisymmetry $f(-u,-v)=-f(u,v)$. For this we may clearly assume that $f\neq 0$. Taking the residues of the Fay identity at $u_1=0$, respectively at $v_1=0$, gives
\begin{align}
a_{-1,0}f(u_2,v_2)+a_{-1,0}f(-u_2,-v_2)&=0,\\
a_{0,-1}f(u_2,v_2)+a_{0,-1}f(-u_2,-v_2)&=0,
\end{align}
and since $(a_{-1,0},a_{0,-1}) \neq (0,0)$, the result follows.
\end{proof}
\begin{rmk} \label{rmk:Fay}
We have already mentioned that our version of the Fay identity is slightly different from the one in \cite[Proposition 5]{BL:MEP}. In particular, the latter does not imply antisymmetry. Indeed, for every $\alpha \neq 0$ the function $f(u,v):=\alpha (\coth(\alpha u)+1)$ satisfies the Fay identity as given in \cite{BL:MEP} (with $u$ corresponding to $\xi$), but does not satisfy Equation \eqref{eqn:Fay}. On the other hand, antisymmetry together with the version of the Fay identity given in \textit{loc.cit.} are equivalent to \eqref{eqn:Fay}. Therefore our version of the Fay identity subsumes antisymmetry as well which is the reason why we prefer to work with it.
\end{rmk}
By Proposition \ref{prop:Faybasic}, we already know that $f=0$, if $(a_{-1,0},a_{0,-1})=(0,0)$. Therefore, the next proposition finishes the proof of Theorem \ref{thm:main} in the special case $a_{-1,0}a_{0,-1}=0$.
\begin{prop} \label{prop:Faybasic2}
Assume that $f \neq 0$ satisfies the Fay identity.
\begin{enumerate}
\item[(i)]If $a_{0,-1}=0$, then $f$ is uniquely determined by $a_{-1,0}$ and $a_{0,1}$, and we have
\begin{equation}
f(u,v)=\frac{a_{-1,0} e^{\gamma uv}}{u}, \quad \mbox{with }\gamma=\frac{a_{0,1}}{a_{-1,0}}.
\end{equation}
\item[(ii)]
If $a_{-1,0}=0$, then $f$ is uniquely determined by $a_{0,-1}$ and $a_{1,0}$, and we have
\begin{equation}
f(u,v)=\frac{a_{0,-1} e^{\gamma uv}}{v}, \quad \mbox{with }\gamma=\frac{a_{1,0}}{a_{0,-1}}.
\end{equation}
\end{enumerate}
\end{prop}
Note that $\gamma$ is well-defined in both cases. Indeed, since $f\neq 0$ we must have $a_{-1,0} \neq 0$ in the first case and $a_{0,-1} \neq 0$ in the second one, by Proposition \ref{prop:Faybasic}.
\begin{proof}
We only prove (i), the proof of (ii) is analogous. Since $a_{0,-1}=0$, Proposition \ref{prop:Faybasic} implies that $f(u,v)$ does not have a pole at $v=0$ and the Fay identity with $v:=v_1=v_2$ and $u:=u_1=u_2$ yields
\begin{equation} \label{eqn:prop2}
f(u,v)^2-2f(u,0)f(2u,v)=0.
\end{equation}
Writing out the Laurent expansion of $f$ and using that $f(u,0)=a_{-1,0}u^{-1}$ (see the proof of Proposition \ref{prop:Faybasic}), Equation \eqref{eqn:prop2} is equivalent to
\begin{equation} \label{eqn:prop3}
\left(\sum_{m,n\geq 0}a_{m,n}u^mv^n\right)^2-a_{-1,0}\sum_{m\geq 1, \, n\geq 0}(2^{m+1}-2)a_{m,n}u^{m-1}v^n=0.
\end{equation}
This implies $a_{m,0}=0$ for $m\neq -1$, and then by induction on $n$ we get $a_{m,n}=0$, if $m \neq n-1$. More generally,
\begin{equation}
\left(\sum_{\genfrac{}{}{0pt}{}{1\leq i,j\leq n-1}{i+j=n}}a_{i-1,i}a_{j-1,j}\right)u^{n-2}v^n-a_{-1,0}\sum_{m\geq 1}(2^{m+1}-2)a_{m,n}u^{m-1}v^n=0,
\end{equation}
for every $n$, showing that $a_{n-1,n}$ is recursively determined by $a_{-1,0}$ and $a_{0,1}$, and therefore $f$ itself is uniquely determined by $a_{-1,0}$ and $a_{0,1}$. On the other hand, by Corollary \ref{cor:solution} there exists a solution to the Fay identity for any given values of $a_{-1,0} \in \bC^{\times}$ and $a_{0,1} \in \bC$ namely
$
\alpha e^{\gamma uv}u^{-1}
$
with $\alpha=a_{-1,0}$ and $\gamma=a_{0,1}/a_{-1,0}$, and this ends the proof.
\end{proof}

\subsection{The ideal of Fay relations}
To study the Fay identity in more detail, it will be convenient to replace the coefficients $a_{m,n} \in \bC$ by symbols $A_{m,n}$ and accordingly to study a formal version of the Fay identity. More precisely, consider the polynomial $\bC$-algebra
\begin{equation} \label{eqn:defA}
\cA:=\bC[\{A_{-1,0},A_{0,-1}\} \cup \{A_{m,n} \, \vert \, m,n \geq 0, \, m+n \mbox{ odd}\}]
\end{equation}
(the restriction on the indices $(m,n)$ is justified by Proposition \ref{prop:Faybasic}) and let
\begin{equation}
\Phi(u,v)=\sum A_{m,n}u^mv^n \in (uv)^{-1}\cA[[u,v]]
\end{equation}
be the generic element where the sum is over all $(m,n)$ as in \eqref{eqn:defA}. By definition it satisfies $\Phi(-u,-v)=-\Phi(u,v)$. Also, define
\begin{equation} \label{eqn:Faygeneric}
\begin{aligned}
\cF(u_1,u_2,v_1,v_2):=\Phi(u_1,v_1)\Phi(u_2,v_2)&+\Phi(-u_2,v_1-v_2)\Phi(u_1+u_2,v_1)\\
&+\Phi(-u_1-u_2,-v_2)\Phi(u_1,v_1-v_2).
\end{aligned}
\end{equation}
A priori, it is contained in $(u_1u_2v_1v_2(u_1+u_2)(v_1-v_2))^{-1}\cA[[u_1,u_2,v_1,v_2]]$.
\begin{lem}
The element $\cF(u_1,u_2,v_1,v_2)$ is contained in $\cA[[u_1,u_2,v_1,v_2]]$, i.e.
\begin{equation}
\cF(u_1,u_2,v_1,v_2)=\sum_{m_1,m_2,n_1,n_2 \geq 0}c_{m_1,m_2,n_1,n_2}u_1^{m_1}u_2^{m_2}v_1^{n_1}v_2^{n_2},
\end{equation}
for some $c_{m_1,m_2,n_1,n_2} \in \cA$.
\end{lem}
\begin{prf}
By definition, $\Phi(u,v)$ has simple poles exactly along $u=0$ and $v=0$ with residues $A_{-1,0}$ and $A_{0,-1}$ respectively. It is therefore sufficient to check that all residues of $\cF$ along $u_1=0$, $u_2=0$, $u_1+u_2=0$, $v_1=0$, $v_2=0$ and $v_1-v_2=0$ vanish, which is straightforward and essentially only uses that $\Phi(-u,-v)=-\Phi(u,v)$.
\end{prf}

\begin{dfn}
Define $\cJ \subset \cA$ to be the ideal generated by the coefficients $c_{m_1,m_2,n_1,n_2}$ of $\cF$. Also, define $\bar\cA:=\cA/\cJ$ to be the corresponding quotient.
\end{dfn}
The ideal $\cJ$ could be called the ideal of Fay relations. By definition, giving a formal Laurent series $f(u,v)=\sum_{m,n>\!\!>-\infty}a_{m,n}u^mv^n \in \bC((u,v))$ which satisfies the Fay identity is equivalent to giving a homomorphism of $\bC$-algebras
\begin{equation}
\varphi_f: \bar\cA \rightarrow \bC, \quad \bar A_{m,n} \mapsto a_{m,n}.
\end{equation}
Solutions to the Fay identity with $a_{-1,0}a_{0,-1} \neq 0$ correspond under this identification to homomorphisms
$
\varphi_f: \bar \cA_0 \rightarrow \bC,
$
where $\bar \cA_0:=\bar \cA \otimes_{\bC} \bC[\bar A^{-1}_{-1,0},\bar A^{-1}_{0,-1}]$. By Proposition \ref{prop:Faybasic2} it is enough to classify the latter, and we shall therefore be interested in understanding the structure of $\bar\cA_0$. To this end, let
$
\bar\cA' \subset \bar\cA$,
be the $\bC$-subalgebra generated by the set $\{\bar A_{0,-1}\} \cup \{\bar A_{m,0} \, \vert \, m \in \{(-1,1,3,5)\}\}$ and denote by
$
\iota: \bar\cA' \hookrightarrow \bar\cA$
the canonical inclusion. Extending scalars to $\bC[\bar A^{-1}_{-1,0},\bar A^{-1}_{0,-1}]$, we get an induced map
\begin{equation}
\iota_0: \bar\cA'_0\rightarrow \bar\cA_0, \quad \mbox{where } \bar\cA'_0:=\bar \cA' \otimes_{\bC} \bC[\bar A^{-1}_{-1,0},\bar A^{-1}_{0,-1}],
\end{equation}
which is clearly injective.
\begin{thm} \label{thm:generators}
The map $\iota_0$ is an isomorphism of algebras.
\end{thm}
More concretely, using the $1$-$1$ correspondence $f \leftrightarrow \varphi_f$ described above, Theorem \ref{thm:generators} says that every solution $f(u,v) \in \bC((u,v))$ to the Fay identity with $a_{-1,0}a_{0,-1} \neq 0$ is uniquely determined by its coefficients $a_{0,-1}$ and $a_{m,0}$, for $m=-1,1,3,5$.
\begin{rmk}
Although we will not need this, one can show that $\bar \cA_0$ is freely generated, as a $\bC[\bar A^{-1}_{-1,0},\bar A^{-1}_{0,-1}]$-algebra, by $\bar A_{0,-1}$ and $\bar A_{m,0}$ for $m=-1,1,3,5$; this follows from Theorem \ref{thm:generators} since for every quintuple $(z_{-1},w_{-1},z_1,z_3,z_5) \in \bC^5$ such that $z_{-1}w_{-1} \neq 0$ there exists a solution $f(u,v)=\sum_{m,n >\!\!> -\infty}a_{m,n}u^mv^n$ to the Fay identity such that $a_{0,-1}=w_{-1}$ and $a_{m,0}=z_m$ for $m=-1,1,3,5$ (see the proof of Theorem \ref{thm:main}).
\end{rmk}

Before we prove Theorem \ref{thm:generators}, we need to introduce some more notation. Let $\fp \subset \cA$ be the ideal generated by all $A_{m,n}$ with $m,n\geq 0$. This is a homogeneous prime ideal of $\cA$, the grading being defined by giving $A_{m,n}$ degree $m+n$. Moreover, since the ideal $\cJ$ of Fay relations is homogeneous, this grading descends to the quotient $\bar \cA$. In general, given a homogeneous ideal $\cI$ of either $\cA$ or $\bar \cA$, we will denote by $\cI_k$ its component of degree $k$.

The following proposition gives explicit formulas for some of the coefficients $c_{m_1,m_2,n_1,n_2}$ and will be the key for proving Theorem \ref{thm:generators}.
\begin{prop} \label{prop:coefficients}
We have the following formulas for the coefficients of $\cF$:
\begin{equation}
c_{0,0,0,0}=3A_{0,-1}A_{0,1}-3A_{-1,0}A_{1,0} \label{eqn:1}
\end{equation}
and for $k\geq 2$ even,
\begin{align}
c_{0,0,0,k}& \equiv (k+3)A_{0,-1}A_{0,k+1}-2A_{-1,0}A_{1,k} \mod \fp^2_k, \label{eqn:2}\\
c_{0,k,0,0}&\equiv 2A_{0,-1}A_{k,1}-(k+3)A_{-1,0}A_{k+1,0} \mod \fp^2_k, \label{eqn:3}\\
c_{0,0,2,k-2}&\equiv \left( \binom{k+2}{3}+1 \right)A_{0,-1}A_{0,k+1}-\left(\binom{k}{2}+\delta_{k,2}\right)A_{-1,0}A_{1,k} \mod \fp^2_k, \label{eqn:4}\\
c_{k-2,2,0,0}&\equiv -\left( \binom{k+2}{3}+1 \right)A_{-1,0}A_{k+1,0}+\left(\binom{k}{2}+\delta_{k,2}\right)A_{0,-1}A_{k,1} \mod \fp^2_k. \label{eqn:5}
\end{align}
Finally, for $0<m<k$ with $k$ as above, we have
\begin{equation}
c_{0,m,0,k-m}\equiv(k+2-m)A_{0,-1}A_{m,k+1-m}-(m+2)A_{-1,0}A_{m+1,k-m} \mod \fp^2_k. \label{eqn:6}
\end{equation}
\end{prop}
\begin{proof}
Since the Fay identity is homogeneous, for every monomial $u_2^rv_2^s$ with $r,s\geq 0$ the coefficient of $u_2^rv_2^s$ in \eqref{eqn:Faygeneric} is congruent modulo $\fp^2_{r+s}$ to the coefficient of $u_2^rv_2^s$ in
\begin{align} \label{eqn:Fayspecial}
\varphi_{r,s}(u_1,v_1)\varphi_{r,s}(u_2,v_2)&+\varphi_{r,s}(-u_2,v_1-v_2)\varphi_{r,s}(u_1+u_2,v_1)\\
&+\varphi_{r,s}(-u_1-u_2,-v_2)\varphi_{r,s}(u_1,v_1-v_2),
\end{align}
where
$
\varphi_{r,s}(u,v)=A_{-1,0}u^{-1}+A_{0,-1}v^{-1}+A_{r,s+1}u^rv^{s+1}+A_{r+1,s}u^{r+1}v^{s}.
$
A straightforward computation of the coefficient of $u_2^mv_2^{k-m}$ in \eqref{eqn:Fayspecial} now yields \eqref{eqn:1}, \eqref{eqn:2}, \eqref{eqn:3} and $\eqref{eqn:6}$.

Similarly, for $k\geq 2$ the coefficient of $v_1^2v_2^{k-2}$ (respectively of $u_1^{k-2}u_2^2$) in \eqref{eqn:Faygeneric} is congruent modulo $\fp^2_k$ to the coefficient of the same monomial in \eqref{eqn:Fayspecial} for $(r,s)=(0,k)$ (respectively for $(r,s)=(k,0)$), and we get \eqref{eqn:4} and \eqref{eqn:5}.
\end{proof}
\begin{proof}[Proof of Theorem \ref{thm:generators}]
It is clearly enough to show that $\bar A_{m,n} \in \bar \cA'_0$ for all $m,n$. We prove this by induction on the degree $d=m+n$ of $\bar A_{m,n}$. For $d=1$, we have $\bar A_{1,0} \in \bar \cA'$ by definition and it follows from \eqref{eqn:1} that $\bar A_{0,1}=\bar A_{-1,0}^{-1}\bar A_{0,-1}\bar A_{1,0} \in \bar \cA'_0$. Note that this also implies that $(\bar\fp^2_0)_2 \subset \bar\cA'_0$ where $\bar\fp^2_0$ denotes the ideal generated by the image of $\fp^2$ in $\bar\cA_0$.

Now we use induction on $d$ to show that $\bar A_{m,n} \in \bar\cA'_0$ for all $m,n$ with $d=m+n$ (this shows in particular that $(\bar\fp^2_0)_{d-1} \subset \bar\cA'_0$). For $d=3$ or $d=5$, since $\bar A_{d,0} \in \bar\cA'$ in that case, we see from \eqref{eqn:3} together with $(\bar\fp^2_0)_{d-1} \subset \bar\cA'_0$ (which follows from the induction hypothesis) that $\bar A_{d-1,1} \in \bar\cA'_0$. Repeating the same argument, using \eqref{eqn:6} and finally \eqref{eqn:2}, we obtain $\bar A_{m,n} \in \bar\cA'_0$ for all $m+n=d$, if $d=3$ or $d=5$.

If $d\geq 7$, the crucial point is that the row vectors $(-(k+3),2)$ and $(-\binom{k+2}{3}-1,\binom{k}{2}+\delta_{k,2})$, where $k=d-1$, are linearly independent. Therefore, from \eqref{eqn:3} and \eqref{eqn:5}, we see that $\bar A_{d,0} \equiv 0 \mod (\bar\fp^2_0)_{d-1}$, i.e. $\bar A_{d,0} \in (\bar\fp^2_0)_{d-1}$. By our induction hypothesis, this shows $\bar A_{d,0} \in \bar\cA'_0$ and using a similar argument as before we get $\bar A_{m,n} \in \bar\cA'_0$ for all $m+n=d$.
\end{proof}

\subsection{Proof of Theorem {\ref{thm:main}}} \label{ssec:3.3}
Let $f(u,v)=\sum_{m,n>\!\!>-\infty}a_{m,n}u^mv^n$ be a solution to the Fay identity. By Proposition \ref{prop:Faybasic}, we have
\begin{equation}
f(u,v)=\frac{a_{-1,0}}{u}+\frac{a_{0,-1}}{v}+\sum_{m,n\geq 0}a_{m,n}u^mv^n, \quad \mbox{with }a_{m,n}=0, \, \mbox{if $m+n \in 2\bZ$.}
\end{equation}
The case $a_{-1,0}a_{0,-1}=0$ has already been taken care of by Propositions \ref{prop:Faybasic} and \ref{prop:Faybasic2} so that we may assume $a_{-1,0}a_{0,-1} \neq 0$, i.e. $a_{-1,0}$ and $a_{0,-1}$ are both invertible. We distinguish between two cases.

\noindent
\textit{Case (i): $(a_{3,0},a_{5,0})=(0,0)$.}
Let
\begin{equation}
\frac{b_{-1,0}}{u}+\frac{b_{0,-1}}{v}+\sum_{m,n\geq 0}b_{m,n}u^mv^n,
\end{equation}
be the Laurent expansion of $e^{\gamma uv}( \alpha/u+\beta/v)$.\footnote{Here and in the following we suppress the dependence of the Laurent coefficients $b_{m,n}$ on the parameters $\alpha,\beta,\gamma,\delta,\tau$.} It is straightforward to verify that setting $\alpha=a_{-1,0}$, $\beta=a_{0,-1}$ and $\gamma=a_{1,0}/\beta$, we have $a_{0,-1}=b_{0,-1}$ and $a_{m,0}=b_{m,0}$ for $m=-1,1,3,5$. Therefore $a_{m,n}=b_{m,n}$ for all $m,n$ by Theorem \ref{thm:generators}. The uniqueness of $\alpha,\beta$ and $\gamma$ is clear.
\bigskip

\noindent
\textit{Case (ii): $(a_{3,0},a_{5,0}) \neq (0,0)$.}
For any choice of parameters $\alpha,\beta,\gamma,\delta,\tau$ as in the statement of Theorem \ref{thm:main}.(ii), we get the following Laurent expansion
\begin{equation}
\delta e^{\gamma uv}F_{\tau}\left(\frac{u}{\alpha},\frac{v}{\beta}\right)=\frac{b_{-1,0}}{u}+\frac{b_{0,-1}}{v}+\sum_{m,n\geq 0}b_{m,n}u^mv^n.
\end{equation}
It follows from \cite[Theorem 3.(iv)]{Zagier:Periods} that $
b_{-1,0}=\alpha\delta$, $b_{0,-1}=\beta\delta,
$ and
\begin{align}
b_{m,0}=\begin{cases}\displaystyle-2\delta\frac{\alpha^{-m}}{m!}G_{m+1}(\tau)=-2b_{-1,0}\frac{\alpha^{-(m+1)}}{m!}G_{m+1}(\tau), & m \neq 1,\\\\
\displaystyle\beta\gamma\delta-2\delta\alpha^{-1}G_2(\tau)=b_{0,-1}\gamma-2b_{-1,0}\alpha^{-2}G_2(\tau), & m=1,
\end{cases}
\end{align}
where $G_k(\tau)$ is the Hecke-normalized Eisenstein series of weight $k$ (the normalization is in fact irrelevant here, we will only need that $G_k(\tau)$ is a modular form of weight $k$).

We now view $b_{3,0}$ and $b_{5,0}$ as functions of $(\alpha,\tau)$ and claim that the map
\begin{equation}
\bC^{\times} \times \overline{\fH} \rightarrow \bC^2 \setminus \{(0,0)\}, \quad (\alpha,\tau) \mapsto (b_{3,0},b_{5,0}),
\end{equation}
where $\overline{\fH}:=\fH \cup \{i\infty\}$, is surjective. Indeed, it is enough to prove that the map $(\alpha,\tau) \mapsto (b_{3,0}^3,b_{5,0}^2)$ is surjective which in turn is equivalent to proving surjectivity of
\begin{equation}
\overline{\fH} \rightarrow \bP^1(\bC), \quad \tau \mapsto [b_{3,0}^3:b_{5,0}^2].
\end{equation}
But this map is surjective because the quotient $b_{3,0}^3/b_{5,0}^2$ is a non-constant modular function (being proportional to $G_4(\tau)^3/G_6(\tau)^2$), and every such function surjects onto $\bP^1(\bC)$.

By what we have just proved, given any coefficients $(a_{3,0},a_{5,0}) \in \bC^2 \setminus \{(0,0)\}$, we can choose $\alpha \in \bC^{\times}$ and $\tau \in \overline{\fH}$ such that we have an equality of Laurent coefficients $b_{3,0}=a_{3,0}$ and $b_{5,0}=a_{5,0}$. Now setting
\begin{equation}
\delta=\frac{a_{-1,0}}{\alpha}, \quad \beta=\frac{a_{0,-1}}{\delta}, \quad \gamma=\frac{a_{1,0}+2a_{-1,0}\alpha^{-2}G_2(\tau)}{a_{0,-1}},
\end{equation}
we get $a_{0,-1}=b_{0,-1}$ as well as $a_{m,0}=b_{m,0}$ for $m=-1,1,3,5$, hence $a_{m,n}=b_{m,n}$ for all $m,n$ by Theorem \ref{thm:generators}.

On the other hand, assume that $\alpha',\beta',\gamma',\delta',\tau'$ are different parameters such that $f(u,v)=\delta'e^{\gamma'uv}F_{\tau'}\left( \frac{u}{\alpha'},\frac{v}{\beta'} \right)$. Then, since the quotient $b_{3,0}^3/b_{5,0}^2$ (which depends neither on $\alpha$, nor on $\alpha'$) is a non-constant modular function, there exists $\left( \begin{smallmatrix}a&b\\c&d\end{smallmatrix} \right) \in \SL_2(\bZ)$ such that $\tau'=\frac{a\tau+b}{c\tau+d}$. The fact that $\alpha'=\alpha(c\tau+d)$ and $\gamma'=\gamma-\frac{c/2\pi i}{\alpha\beta(c\tau+d)}$ then follows by looking at the coefficients $b_{m,0}$ and using the modular transformation property of $G_{m+1}(\tau)$, if $\tau,\tau' \neq i\infty$. Finally, that $\delta'=\frac{\delta}{c\tau+d}$ follows from the previous result using that $\alpha\delta=b_{-1,0}=\alpha'\delta'$ and similarly one shows $\beta'=\beta(c\tau+d)$. This ends the proof of Theorem \ref{thm:main}.
\qed

\section{Proof of the main result} \label{sec:4}

We now give the proof of Theorem \ref{thm:main2}. First, we need to describe the relation between Fay identity and period relations.

\subsection{Comparison with the period relations} \label{ssec:4.1}

We again follow \cite{Zagier:Periods}. Consider the formal generating series
\begin{equation}
C(X,Y;\tau;T)=\frac{(XY-1)(X+Y)}{X^2Y^2}T^{-2}+\sum_{k=2}^{\infty}c_k(X,Y;\tau)\frac{T^{k-2}}{(k-2)!},
\end{equation}
whose coefficients $c_k(X,Y;\tau)$ are given by
\begin{equation}
c_k(X,Y;\tau)=\sum_f \frac{\frac 12(r_f(X)r_f(Y)-r_f(-X)r_f(-Y))}{(2i)^{k-3}\langle f,f\rangle}f(\tau).
\end{equation}
Here the sum is over all normalized Hecke eigenforms $f$ of weight $k$, $r_f(X)$ denotes the extended period polynomial of $f$ and $\langle\cdot,\cdot\rangle$ is the Petersson inner product. It follows from the definition of $C$ and Equation \eqref{eqn:periodrels} that $C$ satisfies the following variant of the period relations
\begin{equation} \label{eqn:periodrelations}
\begin{aligned}
C(X,Y;\tau;T)+C\left(-\frac 1X,Y;\tau;XT\right)&=0,\\
C(X,Y;\tau;T)+C\left(1-\frac 1X,Y;\tau;XT\right)+C\left( \frac{1}{1-X},Y;\tau;(1-X)T \right)&=0.
\end{aligned}
\end{equation}
Using the ``slash operator''
\begin{equation} \label{eqn:action}
(C \vert \gamma)(X,Y;\tau;T):=C\left( \frac{aX+b}{cX+d},Y;\tau;(cX+d)T \right), \quad \gamma=\begin{pmatrix}a&b\\c&d\end{pmatrix} \in \SL_2(\bZ),
\end{equation}
Equation \eqref{eqn:periodrelations} can be expressed more concisely as
$C\vert(1+S)=C\vert(1+U+U^2)=0$, where $S=\left(\begin{smallmatrix}0&-1\\1&0\end{smallmatrix}\right)$ and $U=\left(\begin{smallmatrix}1&-1\\1&0\end{smallmatrix}\right)$ are generators of $\SL_2(\bZ)$ and $\vert$ is extended linearly to the group ring of $\SL_2(\bZ)$.

Now recall from Theorem \ref{thm:Zagier} that $C(X,Y;\tau;T)=F_{\tau}(T,-XYT)F_{\tau}(XT,YT)$. In our context, this gives a relation between the Fay identity and the period relations, which we describe next. Given a formal Laurent series $f(u,v) \in \bC((u,v))$, we define
\begin{equation}
C_f(X,Y,T)=f(T,-XYT)f(XT,YT) \in \bC((X,Y,T)).
\end{equation}
\begin{prop} \label{prop:periodsFay}
The series $C_f$ satisfies the period relations, Equation \eqref{eqn:periodrelations}, if and only if $f$ satisfies the Fay identity, Equation \eqref{eqn:Fay}.
\end{prop}
\begin{proof}
By Proposition \ref{prop:Faybasic}, every solution to the Fay identity is antisymmetric, so we may as well prove that the period relations are equivalent to Fay identity and antisymmetry. The equivalence of antisymmetry for $f$ and the period relation $C_f\vert(1+S)=0$ is straightforward. On the other hand, setting $u_1:=T$, $u_2:=-XT$, $v_1:=-XYT$ and $v_2:=-YT$, the period relation $C_f\vert(1+U+U^2)=0$ is
\begin{equation} \label{eqn:Fayidentity2}
\begin{aligned}
f(u_1,v_1)f(-u_2,-v_2)&+f(-u_2,v_1-v_2)f(-u_1-u_2,-v_1)\\
&+f(-u_1-u_2,-v_2)f(-u_1,v_2-v_1)=0,
\end{aligned}
\end{equation}
which together with antisymmetry implies \eqref{eqn:Fay}. Conversely, Fay identity and antisymmetry together imply \eqref{eqn:Fayidentity2}.
\end{proof}

\subsection{Proof of Theorem {\ref{thm:main2}}} \label{ssec:4.2}

Combining Proposition \ref{prop:periodsFay} and Theorem \ref{thm:main}, we know that $C_f(X,Y,T)$ satisfies the period relations if and only if
\begin{equation}
f(u,v)=e^{\gamma uv}\left(\frac{\alpha}{u}+\frac{\beta}{v}\right), \quad \mbox{ or } \quad f(u,v)=\widetilde{\delta} e^{\gamma uv}F_{\tau}\left( \frac{u}{\alpha},\frac{v}{\beta} \right),
\end{equation}
for some $\alpha,\beta,\gamma \in \bC$ in the first case and for some $\alpha,\beta,\widetilde{\delta} \in \bC^{\times}$, $\gamma \in \bC$, $\tau \in \fH \cup \{i\infty\}$ in the second. With the notation of Theorem \ref{thm:main2}, in the first case we get $C_f(X,Y,T)=C_{P_{\alpha,\beta}}(X,Y,T)$ and in the second case $C_f(X,Y,T)=\delta C(X,Y';\tau;T')$ with $\delta:=\widetilde{\delta}^2$, as claimed. The uniqueness assertion follows from the corresponding statement in Theorem \ref{thm:main}.
\qed


TITLE: l.s.c. on compact subset
QUESTION [1 upvotes]: Let $X$ be a real Banach space. Consider a function $f\colon X\to \mathbb{R}$. Suppose $f$ is l.s.c. . Let $D$ be a compact subset of X. I want to show that $f$ is bounded from below on $D$.
Suppose it's not. That means that there exists a sequance $(x_n)\subset D$ such that $$f(x_{n})\to-\infty.$$ Since $D$ is compact, I can extract a  subsequance $(x_{n_k})$ from $(x_n)$ convergent to some $x_0$. Now I observe that $$-\infty<f(x_0)\le \liminf_{k\to \infty}f(x_{n_k})\le\lim_{k\to \infty}f(x_{n_k})\to-\infty$$
which leads to a contradiction. My question is: how do we know that $\lim_{k\to \infty}f(x_{n_k})\to-\infty$?
My question is motivated by the following argument: we do not know if the sequance $(x_{n})$ is convergent to any point in $D$. We only know that the sequance $f(x_{n})$ tends to $-\infty$. So how do we deduce that $\lim_{k\to \infty}f(x_{n_k})\to-\infty$?

REPLY [1 votes]: Recall that a sequence converges, if and only if, every subsequence converges to the same limit. So if we consider $\{f(x_n)\}$ as a sequence of numbers in $\bar{\mathbb{R}}=[-\infty,\infty]$ we see that by choice $f(x_n)\to -\infty$. In particular, since $\{f(x_{n_k})\}$ is a subsequence of $\{f(x_n)\}$ we have that $\lim_{k\to\infty} f(x_{n_k})=-\infty$ as well.
In your question you mention that we don't know that the sequence $(x_n)$ is convergent to any point in $D$, but this is not important since we are only looking at the values of the $f$ along the subsequence. In particular, the subsequential limit only depends on the limit of the $\mathbb{R}$-valued sequence $\{f(x_n)\}$.


TITLE: Estimating number of 5th power free numbers
QUESTION [1 upvotes]: Let $p_{i}$ be the $i^{th}$ prime. Let $M_{k}$ be the number of $5^{th}$ power free numbers in the set  $\{1, 2, \ldots, (p_{1}^{4}p_{2}^{4}\cdots p_{k}^{4}) \}$, that are divisible by some prime $p \not\in\{p_{1}, p_{2}, \ldots, p_{k}\} $. I am interested in finding $\lim_{k\rightarrow\infty}\frac{M_{k}}{p_{1}^{4}p_{2}^{4}\cdots p_{k}^{4}}$.
I tried using prime number theorem but could not get a useful result. Any help would be greatly appreciated.

REPLY [1 votes]: For any positive integer $q$, we know that the number of integers $n\in\{1,2,\dots,q\}$ that are relatively prime to $q$ equals $\phi(q)$. Since being relatively prime to $q$ depends only on the residue class of $n$ (mod $q$), it follows that for any multiple $Q$ of $q$, the number of integers $n\in\{1,2,\dots,Q\}$ that are relatively prime to $q$ equals $Q\phi(q)/q$. Taking $q=p_1\cdots p_k$ and $Q=q^4$, we see that the proportion of integers $n\in\{1,2,\dots,p_1^4\cdots p_k^4\}$ not divisible by any of $p_1,\dots,p_k$ is exactly $1-\frac{\phi(p_1\cdots p_k)}{p_1\cdots p_k}$, which tends to $0$ as $k$ tends to infinity.
Consequently, $M_k$ is the same as the proportion of integers not divisible by any of $p_{k+1}^5,p_{k+2}^5,\dots$, which one can show approaches $1$ as $k\to\infty$. Indeed, a more careful analysis shows that
$$
M_k = 1 - \frac{\phi(p_1\cdots p_k)}{p_1\cdots p_k} + O(p_{k+1}^{-4}) \sim 1 - \frac1{e^\gamma\log k},
$$
where $\gamma$ is Euler's constant.


TITLE: How to prove $\frac{d}{dx}\sinh x=\cosh x$ when $\sinh$ and $\cosh$ are defined by an integral?
QUESTION [4 upvotes]: Define $\sinh$ and $\cosh$ by
$$x=\int_0^{\sinh x}\frac{dt}{\sqrt{t^2+1}},\, x\in\mathbb{R}$$
$$x=\int_1^{\cosh x}\frac{dt}{\sqrt{t^2-1}},\, x\ge 0$$
and define $\cosh (-x)=\cosh x$ for $x\lt 0$.
By the inverse function theorem, we have
$$\frac{d}{dt}\sinh^{-1}t=\frac{1}{\sqrt{t^2+1}}\implies \frac{d}{dx}\sinh x=\sqrt{\sinh^2 x+1},$$
$$\frac{d}{dt}\cosh^{-1}t=\frac{1}{\sqrt{t^2-1}}\implies \frac{d}{dx}\cosh x=\operatorname{sgn}(x)\sqrt{\cosh^2 x-1}.$$
How can I show that $\frac{d}{dx}\sinh x=\cosh x$ and $\frac{d}{dx}\cosh x=\sinh x$?
I think I should use the identity $\cosh^2 x-\sinh^2 x=1$ but I don't know how to prove that identity from the integral definition. I also managed to prove
$$\frac{d^2}{dx^2}\sinh x=\sinh x$$
but that doesn't seem to really help.

REPLY [3 votes]: $\newcommand{\d}{\,\mathrm{d}}$Notice that if one makes the substitution (for $u\ge0$) $u^2=t^2-1,\,u\d u=t\d t,\,\frac{u}{\sqrt{u^2+1}}\d u=\d t$ then one finds:
$$\begin{align}\int_1^{\sqrt{\sinh^2x+1}}\frac{1}{\sqrt{t^2-1}}\d t&=\int_0^{\sinh x}\frac{u}{u}\frac{1}{\sqrt{u^2+1}}\d u\\&=\int_0^{\sinh x}\frac{1}{\sqrt{t^2+1}}\d t\\&=x\end{align}$$
Which implies that $\cosh x=\sqrt{\sinh^2x+1}=\frac{\d}{\d x}\sinh x$.


TITLE: Determine $\sin(x + y)$ and $\cos(x - y)$, if you know $\cos(x) = \frac{5}{7}, \quad x\in(0,\pi/2), \quad \sin(y) = \frac{1}{5},\quad y\in(\pi/2,\pi)$
QUESTION [1 upvotes]: I need advice on this example.
Determine $\sin (x + y)$ and $\cos (x - y)$, if you know that
$\cos(x) = \frac{5}{7}, \quad x\in(0,\pi/2), \quad \sin(y) = \frac{1}{5}, \quad y\in(\pi/2,\pi).$
I used this formula to express $\sin (x)$ and $\cos (y)$
$$\sin^2(x)+\cos^2(x)=1$$
I recieved:
$$\sin(x)=\sqrt{24}/7\quad\text{and}\quad\cos(y)=\sqrt{24}/5.$$
However, if I now substitute into these formulas I get a different result than I have in the assignment
$$
\begin{align*}
\sin(x+y)&=\sin(x)\cdot \cos(y)+\cos(x)\cdot \sin(y)\\
\cos(x-y)&=\cos(x)\cdot \cos(y)+\sin(x)\cdot \sin(y)
\end{align*}
$$
The correct result should be
$$\cos(x-y) = -8\sqrt{6}/35\quad\text{and}\quad\sin(x+y) = -19/35.$$

REPLY [1 votes]: You have a sign error.
$y$ is in the second quadrant, hence $\cos(y) = -\frac{\sqrt{24}}{5}$.
\begin{align}
\sin(x+y) =  -\frac{\sqrt{24}}{7} \cdot \frac{\sqrt{24}}{5}  + \frac57 \cdot  \frac15 = -\frac{19}{35}
\end{align}


TITLE: Is algorithmic information theory still evolving?
QUESTION [11 upvotes]: I am currently looking for a subject for a thesis and encountered the field of algorithmic information theory.
The field seems very interesting for me, but it seems everything is the field was done before many years.
So my question is:
Is the field “alive” or is it pretty much closed?
Does it have open questions?
Thanks

REPLY [7 votes]: A modern tweak on algorithmic information theory is algorithmic randomness which was developed intensively in the 2000s (2009-2009) and is still quite active.
The most notorious open problem there may be whether Kolmogorov-Loveland randomness (in which martingales are computable but are allowed to bet on bits out of order) is the same as Martin-Löf randomness (in which martingales are only semicomputable, i.e., the capital function is computably approximable from below). This is known to be almost true, e.g. if $A\oplus B=\{2n:n\in A\}\cup \{2n+1:n\in B\}$ is KL-random then either $A$ or $B$ is ML-random.
An example of a recent paper in this area:
Bienvenu, Laurent, Kolmogorov-Loveland stochasticity and Kolmogorov complexity, Theory Comput. Syst. 46, No. 3, 598-617 (2010). ZBL1204.68110..


TITLE: Find the conditional probability density function $f_{y|x}(y,x)$
QUESTION [1 upvotes]: Assuming $Z$ random variable with $Z=X+Y$, $Z$ depending on $X$ but is independent of $Y$.
We know the value of $X$.
Also assume we know the joint pdf $f_{x,z}(x,z)$. Find $f_{y|x}(y,x)$.
Can you give me a hint or solution?
I think we should apply Bayes
$$f_{Y|X}(x,y)=\frac{f_{Y,X}(y,x)}{f_Y(y)}f_{X,Z}(x,z)=f_{X,X+Y}(x,x+y)f_X(x).$$

REPLY [1 votes]: The conditional distribution of $Y$ given $X$ is
$$f_{Y|X}(x,y)=\frac{f_{X,Y}(x,y)}{f_X(x)}.$$
You should be able to find the marginal distribution of $X$ by integrating the joint distribution of $X$ and $Z$:
$$f_X(x)=\int_{-\infty}^\infty f_{X,Z}(x,z)dz.$$
Also, as noted below in the comment by Graham, $f_{X,Y}(x,y)=f_{X,Z}(x,x+y)$. So
$$f_{Y|X}(x,y)=\frac{f_{X,Y}(x,y)}{f_X(x)}=\frac{f_{X,Z}(x,x+y)}{f_X(x)}.$$
In general, consider bivariate functions of random variables: $Z_1=h_1(X,Y)$ and $Z_2=h_2(X,Y)$, and assume we know the joint p.d.f. of $X$ and $Y$, $f_{X,Y}(x,y)$. As long as we can find an inverse transform $X=u(Z_1,Z_2)$ and $Y=v(Z_1,Z_2)$, then the joint p.d.f. of $Z_1$ and $Z_2$ is
$$f_{Z_1,Z_2}(z_1,z_2)=f_{X,Y}\left(u(z_1,z_2),v(z_1,z_2)\right)
\left|
\begin{array}{cc}
\frac{du}{dz_1} & \frac{du}{dz_2}\\
\frac{dv}{dz_1} & \frac{dv}{dz_2}
\end{array}
\right|.$$
For your example: $Z_1=X$ and $Z_2=X+Y$, thus $u(z_1,z_2)=z_1$ and $v(z_1,z_2)=z_2-z_1$. The Jacobian matrix determinant needed is
$$
\left|
\begin{array}{cc}
1 & 0\\
-1 & 1
\end{array}
\right|=1.
$$
So $f_{X,Z}(x,z)=f_{X,Z}(x,x+y)=f_{X,Y}\left(x,y\right)$.


text:
\begin{document}

\title{Advanced Scenario Creation Strategies for Stochastic Economic Dispatch with Renewables}

\author{Ryan~N.~King, Matthew~Reynolds, Devon Sigler, and~Wesley~Jones
\thanks{R. N. King, M. Reynolds, D. Sigler, and W. Jones are with the National Renewable Energy Laboratory (NREL), Golden, CO, 80401 USA, e-mail: ({ryan.king, matthew.reynolds, devon.sigler, wesley.jones}@nrel.gov).}
\thanks{Manuscript received June 6, 2018.}}

\markboth{IEEE Transactions on Power Systems, 2018}
{Shell \MakeLowercase{\textit{et al.}}: Bare Demo of IEEEtran.cls for IEEE Journals}

\maketitle

\begin{abstract}
Real-time dispatch practices for operating the electric grid in an economic and reliable manner are evolving to accommodate higher levels of renewable energy generation.  In particular, stochastic optimization is receiving increased attention as a technique for handling the inherent uncertainty in wind and solar generation.  The typical two-stage stochastic optimization formulation relies on a sample average approximation with scenarios representing errors in forecasting renewable energy ramp events.  Standard Monte Carlo sampling approaches can result in prohibitively high-dimensional systems for optimization, as well as a poor representation of extreme events that challenge grid reliability.  We propose two alternative scenario creation strategies, importance sampling and Bayesian quadrature, that can reduce the estimator's variance.  Their performance is assessed on a week's worth of 5 minute stochastic economic dispatch decisions for realistic wind and electrical system data.  Both strategies yield more economic solutions and improved reliability compared to Monte Carlo sampling, with Bayesian quadrature being less computationally intensive than importance sampling and more economic when considering at least 20 scenarios.
\end{abstract}

\begin{IEEEkeywords}
Stochastic programming, economic dispatch, importance sampling, Bayesian quadrature
\end{IEEEkeywords}

\section*{Nomenclature}
\addcontentsline{toc}{section}{Nomenclature}
\begin{IEEEdescription}[\IEEEusemathlabelsep\IEEEsetlabelwidth{$V_1,V_2,V_3$}]
\item[$x_{g}$]Dispatch of thermal generators
\item[$c_g$]Cost of thermal generator $g$
\item[$L\left(\cdot\right)$]Loss function
\item[$\omega_{w}$]Dispatch of wind generator
\item[$c_w$]Cost of wind plant $w$
\item[$\omega_{w}^{fcst}$]Forecasted wind generation
\item[$\xi_{w}$]Error in wind forecast
\item[$x_{g}^{max}$]Max output of generator $g$
\item[$x_{g}^{min}$]Min output of generator $g$
\item[$R_g^{up}$]Max ramp up of generator $g$
\item[$R_g^{down}$]Max ramp down of generator $g$
\item[$G$]Set of thermal generators
\item[$T$]Set of timesteps
\item[$W$]Set of wind plants
\item[$S$]Set of wind forecast error scenarios
\item[$D$]Set of load buses
\item[$d_{q}$]Value of loads
\item[$y^+_{q}$]Loss-of-load
\item[$y^-_{q}$]Excess capacity
\item[$c^+_{q}$]Cost of Loss-of-load
\item[$c^-_{q}$]Cost of excess capacity
\item[$\omega^{spl}_{w}$]Spilled wind
\item[$c^{spl}_{w}$]Cost of spilled wind
\end{IEEEdescription}

\section{Introduction}
\IEEEPARstart{M}{odern} practices for economic dispatch of the electrical grid are rooted in deterministic problem formulations.  The increasing deployment of stochastic renewable energy generators, such as wind and solar, requires dispatch techniques that are consistent with the reliability and economic operating requirements of the electric grid.  Stochastic optimization techniques have received increased attention from the power systems community, and are a promising way to address the randomness in renewable energy generation.  This paper explores two novel approaches to scenario-based stochastic optimization that yield lower costs and reduced reserve requirements, facilitating the grid integration and further deployment of renewable energy generators.

As access to high performance computing resources grows, grid operators are able to consider increasingly sophisticated optimization problems in their dispatch decisions.  These complexities can include considering AC optimal power flow (ACOPF), larger grids with thousands of buses, multiple stages or time periods, and stochastic optimization techniques.  These techniques can improve grid reliability in the face of uncertain generation, and help remove obstacles to further renewable energy deployment.

An emerging approach to dispatching generators in the face of uncertain renewable generation is to formulate and solve a two-stage stochastic programming problem with recourse.  In two-stage recourse problems a first stage decision is made, followed by a recourse stage where adjustments can be made after a random variable is realized.  In the context of economic dispatch for the grid, this involves making initial generator dispatch decisions that minimize the combined costs of the first stage and the expectation of the recourse costs.  The expected recourse costs are calculated with respect to the distribution of the uncertain renewable generation, and are typically estimated with a sample average approximation using Monte Carlo techniques \cite{constantinescu2011computational,soroudi2013decision}.  However, Monte Carlo techniques exhibit slow convergence with respect to the number of samples and draw samples that are not necessarily a good test of grid robustness in the event of large ramps in generation, creating a need for improved sampling techniques.

In this article we report on two techniques for creating scenarios of errors in forecasting wind generation that improve the convergence rates and accuracy when calculating the expected recourse costs, ultimately achieving more economic dispatch solutions.  The created scenarios also contain ramp events that are more challenging for the grid, resulting in dispatch decisions that produce more reliable grid operation.  The cost savings and reliability improvements are demonstrated in a week-long simulation of 5 minute dispatch decisions on realistic wind and grid data.

In Section II we provide a background on stochastic programming techniques for two-stage recourse problems, the standard sample average approximation, and power grid applications in economic dispatch.  In Section III we introduce new scenario creation strategies using importance sampling and Bayesian quadrature.  In Section IV we demonstrate the implementation of these techniques in a stochastic economic dispatch framework and examine their performance.  Finally, we conclude in Section V with a discussion of the results and extensions to more complex multi-stage problems and more complex network topologies.

\section{Two-Stage Stochastic Programming}
\subsection{General Problem Formulation}
The general two-stage stochastic programming problem \cite{shapiro2009lectures,birge2011introduction} is defined as:
\begin{equation}
\min_{\bm{x}} \quad f\left(\bm{x}\right) + \mathbb{E}_{\bm{\xi}}\left[Q\left(\bm{x},\bm{\xi}\right)\right]\,,
\end{equation}
where $\bm{x}$ are the first stage decision variables, $Q\left(\bm{x},\bm{\xi}\right)$ is the second-stage cost function defined by,
\begin{align}
Q\left(\bm{x,\xi}\right) =  \min_{\bm{y}} \quad & q\left(\bm{x},\bm{\xi}\right)  \nonumber \\
\text{s. t.} \quad & \bm{T_{\xi}x + W_{\xi}y=b_{\xi}},
\end{align}
and $\bm{y}$ are the second-stage variables, and $\bm{\xi}$ are random variables.  General constraint equations are represented by $\bm{T_{\xi}}$, $\bm{W_{\xi}}$ and $\bm{b_{\xi}}$, and will be made specific to grid applications in the next section.  We restrict ourselves to using the linear stochastic programming problem, i.e.
\begin{equation}
\min_{\bm{x}} \quad \bm{c}^T\bm{x} + \mathbb{E}_{\bm{\xi}}\left[L\left(\bm{x},\bm{\xi}\right)\right]\,,
\end{equation}
where the recourse function is defined as the solution to,
\begin{align}
L\left(\bm{x,\xi}\right) =  \min_{\bm{y}} \quad & \bm{c_{\xi}}^{\sf T}\bm{y}  \nonumber \\
\text{s. t.} \quad & \bm{T_{\xi}x + W_{\xi}y=b_{\xi}}\\
& \bm{y}\geq \bm{0} \,. \nonumber
\end{align}
In the context of economic dispatch, the first stage variables $\bm{x}$ represent the dispatch levels of all generators, $\bm{c}$ are linear costs, and the expectation $\mathbb{E}_{\bm{\xi}}\left[\cdot\right]$ is taken with respect to the random variable $\bm{\xi}$, which represents errors in forecasting the available wind generation.

The second stage recourse variable, $\bm{y}$, represents adjustments to generation using reserves, such as up/down regulation, flex reserve or automatic generation control (AGC), that are necessary based on the amount of wind generation that materializes.

The optimization problem embedded inside the expectation $\mathbb{E}_{\bm{\xi}}\left[\cdot\right]$ makes it difficult to solve, so a sample average approximation (SAA) is commonly used to represent the expectation.  In the SAA approach, we introduce a set of $N$ scenarios, $\mathcal{S} = \left\lbrace\bm{\xi}_i\right\rbrace_{i=1}^N$, and approximate the expectation as $\mathbb{E}_{\bm{\xi}}\left[L\left(\bm{x},\bm{\xi}\right)\right] \approx \frac{1}{N}\sum_{i=1}^{N}L\left(\bm{x},\bm{\xi}_i\right)$.  Commuting the summation and minimization, we can simplify to the extensive form of a two-stage stochastic programming problem:
\begin{eqnarray}
\min_{\bm{x,y}_i i=1,\ldots,N} \quad \bm{c}^{\sf T}\bm{x} + \frac{1}{N}\sum_{i=1}^{N}\bm{c}_{\bm{\xi}_i}^{\sf T}\bm{y}_i \hspace{1.55in} & \nonumber \\
\begin{aligned}
\text{s. t.} \quad & \bm{Ax} & & & & & &=\bm{b} \\
& \bm{T}_{\bm{\xi}_1}\bm{x} \quad + & \bm{W}_{\bm{\xi}_1}\bm{y}_1 & & & & &= \bm{b}_1 \\
& \bm{T}_{\bm{\xi}_2}\bm{x} \quad + & & \bm{W}_{\bm{\xi}_2}\bm{y}_2  & & & &= \bm{b}_2 \\
& \bm{T}_{\bm{\xi}_3}\bm{x} \quad + & & & \ddots  & & &= \vdots \\
& \bm{T}_{\bm{\xi}_N}\bm{x} \quad + & & &  & \bm{W}_{\bm{\xi}_N}\bm{y}_N & &= \bm{b}_N
\\
& \bm{x}\geq\bm{0}, & \bm{y}_1 \geq \bm{0}, & \quad \bm{y}_2 \geq \bm{0}, & \ldots \:, & \bm{y}_N\geq\bm{0} \,.
\end{aligned}
\end{eqnarray}

A major challenge in the SAA approach is selecting the scenarios $\mathcal{S}$ such that the expectation $\mathbb{E}_{\bm{\xi}}\left[L\left(\bm{x},\bm{\xi}\right)\right]$ is accurately and efficiently calculated, while also ensuring reliable operation of the grid during different levels of wind generation.  In Monte Carlo approaches, the samples are drawn from their nominal distribution $p\left(\bm{\xi}\right)$ and given equal weights.  This approach yields $\mathcal{O}\left(N^{-1/2}\right)$ convergence, requiring prohibitively many samples to achieve an acceptable error level.  Other methods for selecting scenarios have been recently proposed, including scenarios motivated by grid reserve requirements \cite{papavasiliou2011reserve}, chance constraints \cite{cheng2018chance,lubin2016robust}, polynomial chaos expansions \cite{safta2017efficient}, or epi-spline basis functions \cite{staid2017generating}.

\subsection{Two-Stage Stochastic Economic Dispatch}
In this is study we restrict ourselves to studying the stochastic economic dispatch problem, i.e. balancing the (stochastic) demand for power across a grid with the (stochastic) power generation from all generators on the network. This study considers aggregated wind and load data, and therefore our model is formulated without network constraints. However, the scenario creation methods presented in this study should generalize to more complicated constraints on the network, e.g. DC and ACOPF. For a review of the economic dispatch problem, as well as the DCOPF and ACOPF problems, we recommend \cite{wood2012power,zhu2015optimization,taylor2015convex}. For an introduction to stochastic grid operations we recommend \cite{conejo2010decision}.

In the context of grid operations, the stochastic economic dispatch problem takes the following form:
\begin{align}
\min_{\substack{\bm{x} }} \quad & \sum_{\substack{g\in G}} \left( c_g x_{g} +\mathbb{E}_{\bm{\xi}}\left[L\left(\bm{x},\bm{\xi}\right)\right] \right) \nonumber  \\
\text{s. t.} \quad & x_g^{min} \leq x_{g} \leq x_g^{max} \quad \forall \: g \\
& -R_g^{down}\leq x_g - I_g \leq R_g^{up} \quad \forall \: g, \nonumber
\end{align}
where $\bm{x} = x_{g}$ for $g$ in the set of generators, $G$, are the first stage decisions representing thermal dispatch, and $L\left(\bm{x},\bm{\xi}\right)$ is the recourse function for a particular realization of the random variable $\bm{\xi}$, which represents the error in forecasting wind generation at each wind plant $w$ in the set of wind plants $W$.  The recourse function is defined as the solution to the following optimization problem:
\begin{align}\label{eq:loss-function}
&L\left(\bm{x,\xi}\right) =  \nonumber \\
\min_{\substack{\bm{y}^{\pm},\bm{\omega},\bm{\omega}^{spl} }}\quad  &\sum_{\substack{w\in W}} \left(c_w \omega_{w} + c_w^{spl} \omega_{w}^{spl} \right) \quad + \nonumber \\
& \sum_{g\in G} \left( c_q^+ y^+_{q} + c_q^- y^-_{q} \right) \quad + \nonumber \\
\text{s. t.} \quad
& 0\leq y^+_{q} \quad \forall \: q \\
& 0\leq y^-_{q} \quad \forall \: q \nonumber \\
& 0  \leq  \omega_{w}  \leq \omega^{fcst}_{w} + \xi_{w}  \quad \forall \: w \nonumber \\
& \omega_{w}^{spl} = \left(\omega^{fcst}_{w} + \xi_{w}\right) \nonumber \\
& \qquad - \left( \omega_{w} \right)  \quad \forall \: w \nonumber \\
& \sum_{q \in D} \left(y^+_{q} - y^-_{q}\right) + \sum_{w \in W} \omega_{w} + \sum_{g\in G} x_{g}  \nonumber \\
& \qquad =\sum_{q \in D}\left(d_{q}  \right), \nonumber
\end{align}
where $\bm{\xi}=\xi_w, \, \bm{\omega}=\omega_w, \, \bm{\omega}^{spl}=\omega^{spl}_w \text{ for } w\in W$ and $\bm{y}^{\pm}=y^{\pm}_q$ for each bus $q$ in the set of buses $D$. The expectation $\mathbb{E}_{\bm{\xi}}\left[\cdot\right]$ is taken with respect to the random variable $\bm{\xi}$.  Generator ramp constraints are imposed between sequential timesteps, and not allowed to participate in recourse.  The loss of load or excess capacity represented by $\bm{y}^{\pm}$ is a slack variable that facilitates finding feasible solutions and represents a challenge to grid reliability and incurs substantial costs.

We assume the wind forecast errors are relative to a persistence forecast based on how well this approach performs in forecasting over short time-periods (see e.g. \cite{giebel2003state}). For an in-depth survey of wind power forecasting methods we recommend \cite{foley2012current} to the reader.

\section{Strategies for Scenario Generation}

\subsection{Importance Sampling}
Solving the two-stage stochastic economic dispatch problem requires a method of selecting scenarios. A logical starting point would be to draw from the nominal distribution of forecast errors in the wind power. However, this approach will not take into account costs communicated by the loss function \eqref{eq:loss-function}, potentially resulting in scenarios that do not include information about expensive events. As an example, consider a single wind farm generating power for one time period. Assuming that the nominal distribution of forecast errors is symmetric about zero, simply sampling from this distribution will produce scenarios that represent positive forecast errors, i.e. more wind than forecast, about half the time. If we use a small number of scenarios, e.g. $5$ to $10$, as is common in stochastic optimization problems, the scenario set may not include scenarios that represent negative forecast errors. This would be a costly mistake since both loss-of-load and reserves are more expensive than dispatching extra generation during the first stage decision process.

To address this issue with random sampling we employ importance sampling, see e.g. \cite{bucklew_introduction_2004}. The reason for using importance sampling is that by including information from the cost function $L\left(\bm{x},\bm{\xi}\right)$, we construct a new probability density function that assign more mass to values of $\bm{\xi}$ corresponding to errors in overpredicting the available wind generation. Additionally, the construction of the importance distribution also provides a formula for properly weighting these samples to ensure the weighted average actually converges to $\mathbb{E}_{\bm{\xi}}\left[L\left(\bm{x},\bm{\xi}\right)\right]$.  Finally, importance sampling can also yield a lower variance estimate of the expected recourse cost.

Recall that our expected recourse cost is
\begin{equation}
\mu = \mathbb{E}_{\bm{\xi}}\left[L\left(\bm{x},\bm{\xi}\right)\right]=\int_{\Omega} L\left(\bm{x},\bm{\xi}\right)p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi}
\end{equation}
where $p$ is a probability density function on $\mathbb{R}^d$ and $p\left(\bm{\xi}\right)=0$ for all $\bm{\xi}\notin \Omega$.  Introducing a new pdf $q$ on $\mathbb{R}^d$ (the importance distribution) we can write
\begin{align}
\mu = \mathbb{E}_{\bm{\xi}}\left[L\left(\bm{x},\bm{\xi}\right)\right] & =\int_{\Omega} L\left(\bm{x},\bm{\xi}\right)p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi} \\
& = \int_{\Omega} \frac{L\left(\bm{x},\bm{\xi}\right)p\left(\bm{\xi}\right)}{q\left(\bm{\xi}\right)}q\left(\bm{\xi}\right)\mathrm{d}\bm{\xi} \\
& = \mathbb{E}_q\left[\frac{L\left(\bm{x},\bm{\xi}\right)p\left(\bm{\xi}\right)}{q\left(\bm{\xi}\right)}\right]
\end{align}
where $\mathbb{E}_q\left[\cdot\right]$ denotes expectation with $\bm{\xi}$ distributed as $q$.  The ratio $p\left(\bm{\xi}\right)/q\left(\bm{\xi}\right)$ is an adjusted weight on sample evaluations of $L$ that reflects sampling from $q$ instead of $p$.

The importance sampling estimate of $\mu$ is then
\begin{equation}
\widehat{\mu}_q = \frac{1}{N}\sum_{i=1}^{N}\frac{L\left(\bm{x},\bm{\xi}_i\right)p\left(\bm{\xi}_i\right)}{q\left(\bm{\xi}_i\right)},
\end{equation}
where $\bm{\xi}_i$ are realizations of the random variable $\bm{\xi}$ with distribution $q\left(\bm{\xi}\right)$. The optimal importance sampling distribution, defined as $q\left(\bm{\xi}\right)=L\left(\bm{x},\bm{\xi}\right)p\left(\bm{\xi}\right)/\mu$, is a zero-variance estimator, and computes the correct expectation from a single sample.  Such an estimator suffers from the \textit{curse of circularity} (the quantity being estimated, $\mu$, is required to properly normalize the importance distribution), but provides guidance on constructing an approximate importance sampling distribution.

Our approach to constructing an importance distribution is to solve a deterministic economic dispatch problem to approximate $\mu$ and $L\left(\bm{x},\bm{\xi}\right)$, and then use them to build an approximation of $q\left(\bm{\xi}\right)$. Our first simplification involves approximating the loss function as $L\left(\bm{x},0\right)$, i.e. we assume that the wind forecast is perfect. This step reduces the two stage stochastic optimization problem to a simple linear program:
\begin{align}
\bm{x}^* \: = \: \arg\min_{\substack{\bm{x} }} \quad & \sum_{\substack{g\in G}} \left( c_g x_{g} +L\left(\bm{x},\bm{0}\right) \right) \nonumber  \\
\text{s. t.} \quad & x_g^{min} \leq x_{g} \leq x_g^{max} \quad \forall \: g \\
& -R_g^{down}\leq x_g - I_g \leq R_g^{up} \quad \forall \: g, \nonumber
\end{align}
which we solve for $\bm{x}$. Given the solution $\bm{x}^*$ to the simplified dispatch problem we compute the expectation $\tilde{\mu}=\mathbb{E}\left[L\left(\bm{x}^*,\bm{\xi}\right)\right]$ via the trapezoidal rule. We remark that this process is fast since each evaluation of the loss function $L\left(\bm{x}^*,\bm{\xi}\right)$ is computed by solving a linear program. Given the approximated expected cost, $\tilde{\mu}$, the importance sampling distribution is defined as
\begin{equation}
q\left(\bm{\xi}\right)=\frac{L\left(\bm{x}^*,\bm{\xi}\right)}{\tilde{\mu}}p\left(\bm{\xi}\right)\,,
\end{equation}
and $q\left(\bm{\xi}\right)$ is then sampled instead of $p\left(\bm{\xi}\right)$ in the solution of the original stochastic economic dispatch problem.

\subsection{Gaussian Process Representation}
One shortcoming in the importance sampling approach presented in the preceding sub-section is that in reality, $\bm{x}$ is a first-stage variable in a two-stage stochastic optimization problem, and its value is computed at the same time as the solution to the rest of the problem. By ignoring this coupling, our proxy $\bm{x}^*$ may be off by a significant amount, leading to an importance distribution with an incorrect shape and location.

We propose an alternative to using the importance sampling approach: representing the loss function as a Gaussian process (GP) and computing its expectation via a Bayesian quadrature technique. One benefit of this approach is that handles uncertainty in the second stage loss function introduced by not knowing the first stage dispatch.  In other words, we can treat $L\left(\bm{\xi}\right)$ as a random function when $\bm{x}$ is unknown.  Additionally, computing expectations of the GP with Bayesian quadrature does not require knowledge of the loss function to select the quadrature points. Instead, the Bayesian quadrature algorithm spreads out the sample points $\bm{\xi}_i$ in the domain based on the nominal distribution $p\left(\bm{\xi}\right)$ and the Gaussian process covariance kernel function.  The Bayesian quadrature algorithm also learns quadrature weights, thus ensuring that the weighted sum of the loss function sampled at the quadrature points returns the correct value for the expectation loss. A comprehensive review of Gaussian processes is provided by Rasmussen and Williams \cite{rasmussen_gaussian_2006}, while a description of the Bayesian quadrature sampling strategy can be found in \cite{ghahramani_bayesian_2003}.

By representing the loss function as a GP, we make the assumption that the joint distribution of any finite number of loss function evaluations is a multivariate Gaussian.  The GP is fully specified by a mean function $m\left(\bm{\xi}\right) = \mathbb{E}\left[L\left(\bm{\xi}\right)\right]$ and a kernel function $k\left(\bm{\xi},\bm{\xi}'\right) = \mathbb{E}\left[\left(L\left(\bm{\xi}\right)-m\left(\bm{\xi}\right)\right)\left(L\left(\bm{\xi^\prime}\right)-m\left(\bm{\xi^\prime}\right)\right)\right]$.  When represented as a GP, the loss function is
\begin{equation}
L\left(\bm{\xi}\right) \sim \mathcal{GP} \left[m\left(\bm{\xi}\right),k\left(\bm{\xi},\bm{\xi}'\right)\right]\,,
\end{equation}
and returns a Gaussian distribution when evaluated at a particular value $L\left(\bm{\xi}^*\right)$.
The kernel function describes the covariance between any two points in the GP and encodes structural properties of the underlying random function like the smoothness, periodicity, and stationarity.  The kernel function is a nonlinear measure of similarity between inputs, with the dot product kernel leading to standard linear regression.  In this study we use a common nonlinear kernel function, the Gaussian or squared exponential function:
\begin{equation}
k\left(\bm{\xi},\bm{\xi}'\right) = \tau^2\exp\left(-\frac{|\bm{\xi}-\bm{\xi}'|^2}{2l^2}\right)
\end{equation}
where $\tau$ describes the magnitude and $l$ describes the smoothness of the function or correlation length between two points.

The GP can be conditioned on a set of loss function training points $\bm{\Xi}=\left\lbrace\bm{\xi}_i\right\rbrace_{i=1}^N$ to improve the posterior estimate of the loss function at new test points.  The joint distribution of training outputs $\mathbf{L}=L\left(\bm{\Xi}\right)$ and a test output $\mathbf{L}^*=L\left(\bm{\xi}^*\right)$ is
\begin{equation}
\left[\begin{array}{c} \mathbf{L} \\ \mathbf{L}^*\end{array}\right] \sim
\mathcal{N}\left(\left[\begin{array}{c} m\left(\bm{\Xi}\right) \\ m\left(\bm{\xi}^*\right)\end{array}\right], \left[\begin{array}{cc} k\left(\bm{\Xi},\bm{\Xi}\right) & k\left(\bm{\Xi},\bm{\xi}^*\right) \\ k\left(\bm{\xi}^*,\bm{\Xi}\right) & k\left(\bm{\xi}^*,\bm{\xi}^*\right)\end{array}\right] \right) \,.
\end{equation}
Conditioning on the training data $\mathcal{D}=\left \lbrace\bm{\Xi},\mathbf{L}\right\rbrace$ yields the following estimate for $\mathbf{L}^*$:
\begin{gather}
L\left(\bm{\xi}^* | \mathcal{D}\right)  \sim \mathcal{N}\left[m_{\mathcal{D}}\left(\bm{\xi}^*\right), \text{cov}_{\mathcal{D}}\left(L\left(\bm{\xi}^*\right),L\left(\bm{\xi}\right)\right) \right]\\
m_{\mathcal{D}}\left(\bm{\xi}^*\right) = \mathbf{k}_*^{\sf T}\mathbf{K}^{-1}\mathbf{L} \\
\text{cov}_{\mathcal{D}}\left(L\left(\bm{\xi}^*\right),L\left(\bm{\xi}\right)\right) = k\left(\bm{\xi}^*,\bm{\xi}^*\right) - \mathbf{k}_*^{\sf T}\mathbf{K}^{-1}\mathbf{k}_*
\end{gather}
where we have introduced the shorthand $\mathbf{k}_*=k_{*,i}=k\left(\bm{\Xi}_i,\bm{\xi}^*\right)$ and $\mathbf{K}=K_{ij}=k\left(\bm{\Xi}_i,\bm{\Xi}_j\right)$.  Note that the mean prediction $m_{\mathcal{D}}\left(\bm{\xi}^*\right)$ is simply a linear combination of the training points $\mathbf{L}$ and requires inverting a matrix that scales with the number of training points, N.  Additionally, the posterior covariance is given by the prior covariance less the information gained from the training points.

\subsection{Bayesian Quadrature}
Recall that in stochastic economic dispatch we are faced with calculating the expectation of our GP loss function $\mu = \int_{\Omega} L\left(\bm{\xi}\right)p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi}$.  In Bayesian quadrature \cite{ghahramani_bayesian_2003,ohagan_bayeshermite_1991,hennig_probabilistic_2015} the integral $\mu$ is modeled as a Gaussian random variable because integration is a linear transformation that preserves Gaussianity and the integrand, $L\left(\bm{\xi}\right)$, is a Gaussian process.  The density function $p\left(\bm{\xi}\right)$ must also be a Gaussian, however this can be achieved by re-weighting via importance sampling.

Bayesian quadrature allows us to express the expectation and variance of  the loss analytically, and then devise a sampling strategy to optimally reduce this variance.  Following the derivation in \cite{ghahramani_bayesian_2003}, the loss averaged over possible functions is simply the expectation of the mean GP loss function:
\begin{align}
\mathbb{E}_{L|\mathcal{D}} \left[\mu\right] &= \int\int L\left(\bm{\xi}\right)p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi} \; p\left(L|\mathcal{D}\right) \mathrm{d}L\\
& = \int \left[ \int L\left(\bm{\xi}\right) p\left(L|\mathcal{D}\right) \mathrm{d}L\right] p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi} \\
& = \int m_{\mathcal{D}}\left(\bm{\xi}\right)p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi}\\
& = \int k\left (\bm{\xi},\bm{\Xi}\right)\mathbf{K}^{-1}\mathbf{L} \, p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi}\,.
\end{align}
Similarly, the expected variance of the estimate with respect to possible functions is
\begin{align}
\mathbb{V}_{L|\mathcal{D}} \left[\mu\right] &=
 \int\left[ \int L\left(\bm{\xi}\right)p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi} - \right. \nonumber \\
& \ \ \  \left. \int m_{\mathcal{D}}\left(\bm{\xi}'\right)p\left(\bm{\xi}'\right)\mathrm{d}\bm{\xi}' \right]^2 p\left(L|\mathcal{D}\right) \mathrm{d}L \\
&=
\int\int\int\left[ L\left(\bm{\xi}\right) - m_{\mathcal{D}}\left(\bm{\xi}\right)\right] \left[ L\left(\bm{\xi}'\right) - \right. \nonumber\\
& \ \ \ \left. m_{\mathcal{D}}\left(\bm{\xi}'\right)\right]p\left(L|\mathcal{D}\right) \mathrm{d}L \; p\left(\bm{\xi}\right) p\left(\bm{\xi}'\right)\mathrm{d}\bm{\xi}\mathrm{d}\bm{\xi}'  \\
&=\int \int \text{cov}_{\mathcal{D}}\left(L\left(\bm{\xi}\right),L\left(\bm{\xi}'\right)\right)p\left(\bm{\xi}\right) p\left(\bm{\xi}'\right)\mathrm{d}\bm{\xi}\mathrm{d}\bm{\xi}'\\
&= \int\int \left[ k\left(\bm{\xi},\bm{\xi}'\right) - \right. \nonumber\\
& \ \ \ \left. k\left (\bm{\xi},\bm{\Xi}\right) \, \mathbf{K}^{-1} \, k\left(\bm{\Xi},\bm{\xi}\right) \right]p\left(\bm{\xi}\right) p\left(\bm{\xi}'\right)\mathrm{d}\bm{\xi}\mathrm{d}\bm{\xi}'\,.
\end{align}

A useful sampling strategy is to choose the set of samples that minimizes the expected variance:
\begin{equation}
\bm{\Xi}^* = \text{arg}\min_{\bm{\Xi}} \quad \mathbb{V}_{L|\mathcal{D}} \left[\mu\right] \,.
\end{equation}
This  can be compactly represented as
\begin{equation}
\bm{\Xi}^* = \text{arg}\min_{\bm{\Xi}} \quad Z - \mathbf{w}^{\sf T} \mathbf{K}^{-1}\mathbf{w} \,,
\end{equation}
where we have defined the following quantities
\begin{gather}
\mathbf{w}=w_i := \int k\left(\bm{\xi},\bm{\Xi}_i\right)p\left(\bm{\xi}\right)\mathrm{d}\bm{\xi} \\
 Z = \int \int k\left(\bm{\xi},\bm{\xi}'\right) p\left(\bm{\xi}\right)p\left(\bm{\xi}'\right)\mathrm{d}\bm{\xi}\mathrm{d}\bm{\xi}'\,.
\end{gather}
Because the integrated prior covariance $Z$ has no dependence on the sample points, the optimization can be simplified to
\begin{equation} \label{eq:BQopt}
\bm{\Xi}^* = \text{arg}\max_{\bm{\Xi}} \quad\mathbf{w}^{\sf T} \mathbf{K}^{-1}\mathbf{w} \,,
\end{equation}
which shows that the optimal sampling strategy only depends on the choice of kernel covariance function $k\left(\bm{\xi},\bm{\xi}'\right)$ and the nominal distribution $p\left(\bm{\xi}\right)$, but not on the actual loss function evaluations $L\left(\bm{\Xi}\right)$.  To develop some intuition for Eq. \ref{eq:BQopt}, we can see that maximization involves a balance between selecting points that reflect the nominal distribution $p\left(\bm{\xi}\right)$ while avoiding redundancies from selecting similar points.

That the sampling strategy is independent of the loss function evaluations is a key differentiator between importance sampling and Bayesian quadrature approaches.  This is an advantage for Bayesian quadrature in cases where the loss function is unknown or stochastic, e.g. due to uncertainties in electrical loads.  Furthermore, as a quadrature rather than sampling strategy, the Bayesian quadrature points can be reused at successive timesteps if $k\left(\bm{\xi},\bm{\xi}'\right)$ and $p\left(\bm{\xi}\right)$ are unchanged.  The Bayesian quadrature estimate of the expected loss function takes the form of the following quadrature rule

\begin{equation}
\mathbb{E}_{L|\bm{\Xi},\bm{L}} \left[\mu\right] \approx \mathbf{w}^{\sf T}\mathbf{K}^{-1}\mathbf{L}=\sum^N_{i=1}\gamma_i L\left(\bm{\Xi}_i\right)
\end{equation}

\section{Performance of Different Sampling Strategies}

\subsubsection{Test data}

To study the performance of our economic dispatch methods we examine a flattened version of the Reliability Test System, modified by the Grid Modernization Lab Consortium (RTS-GMLC) \cite{GMLCRTS}. The RTS-GMLC system is an update to the original 1979 \cite{chairman1979ieee} and 1996 \cite{grigg1999ieee} version of the RTS that contains 73 buses, 96 conventional thermal generators and 4 wind plants. However, for our experiments we flatten the network and only consider power balance, i.e. matching the aggregated power output from the 96 conventional thermal generators and 4 wind plants to the aggregated demand without taking into consideration the network

The wind plant data used in our experiments is drawn from the Wind Integration National Dataset (WIND) Toolkit \cite{DRAXL2015355}. To appropriately model wind farms using WIND Toolkit data, wind plant power outputs were aggregated by first assigning geographic locations to nodes in the RTS-GMLC system and then summing the output power from the wind plants nearest to the node locations until the prescribed annualized capacity for the node in the RTS-GMLC system was reached.

The experiments in this section take place during a week in July. Load data from the RTS-GMLC system was combined with simulated wind data from 2013 from the WIND Toolkit. To determine $p\left(\bm{\xi}\right)$, we fit a Student's t-distribution to errors in a 5 minute persistence forecast.  A t-distribution was chosen because of the non-normal, heavy tailed behavior of changes in wind generation. The observation heavy-tailed distributions provide a good fit to wind data has been noted previously, see, e.g. \cite{hodge2011wind}.  The distributions of wind forecast errors could be further refined by conditioning on current power, time of day, season, etc, however this study used a fixed t-distribution that best matched errors aggregated over the full month of July.

The actual loads and wind production data used in our experiments are shown in Figure~\ref{fig:load_and_generation}. There are two important aspects of the data: first, the peak renewable generation is just over one quarter of the peak demand, corresponding to a moderate renewable energy penetration scenario. Second, there is a wind ramp event that occurs around midnight July 25th. The wind ramp event will pose an economic and reliability challenge to test the different sampling strategies in our experiments.

\begin{figure}[!t]
\centering
\includegraphics[width=3.0in]{load_and_gen.png}
\caption{Wind and load actuals used in economic dispatch experiments. We make note of the wind ramp event occurring on July 26th.}
\label{fig:load_and_generation}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[width=3.5in]{sampling.png}
\caption{The upper left panel shows the loss function (red), nominal distribution of $\xi$ (green), and importance distribution (blue) of forecast errors.  The other panels show representative samples generated by each sampling technique.}
\label{fig:importance-dist}
\end{figure}

We examined the performance of each sampling strategy with 5, 10, 20, and 50 scenarios drawn at each 5 minute dispatch decision interval.  The resulting costs of the first and second stage decisions for the full week are shown in Figure ~\ref{fig:scenarioTimeseries} for all three strategies and the four sets of scenarios.  Additionally, a summary of the costs for the full week is shown in Table~\ref{tab:costs}.  To make the differences in sampling strategies more clear, our model does not explicitly implement reserves.  Instead any shortfall in generation results in an expensive loss-of-load penalty incurred during the second stage.  This results in the bulk of the costs attributed to first stage dispatch decisions, and rare, but comparatively large second stage costs.  In practice, these loss-of-load events would largely be addressed by ISO policies on reserve requirements.

Representative samples produced by each sampling technique, along with the nominal distribution, importance distribution, and loss function, are shown in Figure \ref{fig:importance-dist}.   The outcomes of these sampling strategies are shown in Figure \ref{fig:scenarioTimeseries} which displays a week-long timeseries of first and second stage dispatch cost.  The costs broken down by stage are additionally summarized in Table \ref{tab:costs2} and Table \ref{tab:costs3}.

The Monte Carlo method that draws from the nominal distribution has the most loss-of-load events. Alternatively, drawing from the importance sampling distribution leads to the smallest number of loss-of-load events. The reason for this behavior is clear: the shape of the loss function, large penalties for loss-of-load and relatively small penalties for overload or wind spilling, causes the importance distribution to favor sampling scenarios of wind overprediction. Thus, we are more likely to draw scenarios that represent lower than forecast values for wind from the importance distribution. These scenarios allow the stochastic optimization algorithm to hedge against expensive loss events. Drawing from the nominal distribution with Monte Carlo sampling, however, leads to scenarios that represent up and down ramps in wind power generation with equal probability. Therefore, for many of the 5 minute time-periods, the stochastic optimization algorithm was not equipped to hedge against loss-of-load events with Monte Carlo samples.

The Bayesian quadrature approach results in a consistent decrease in maximum second stage costs as the number of scenarios is increased.  This consistent performance improvement is a strength of a quadrature-based approach.  In contrast, importance sampling does not show a monotonic decrease in peak second stage costs as more scenarios are considered.  This is likely due to errors arising from approximating the true importance distribution.

\begin{table}[!t]
\caption{Economic Dispatch Total Costs}
\label{tab:costs}
\centering
\begin{tabular}{cccc}
\hline
& \multicolumn{3}{c}{Sampling Strategy Costs (\$)} \\
\cline{2-4}
\# of Scenarios    & MC & IS & BQ\\
\hline
5      & \num{1.453e7}   & \num{1.383e7} & \num{1.394e7}      \\
10       & \num{1.409e7}   & \num{1.382e7} & \num{1.394e7}      \\
20       & \num{1.391e7}   & \num{1.383e7} & \num{1.378e7}      \\
50      & \num{1.386e7}   & \num{1.382e7} & \num{1.378e7}      \\
\hline
\end{tabular}
\end{table}

\begin{table}[!t]
\caption{Economic Dispatch First Stage Costs}
\label{tab:costs2}
\centering
\begin{tabular}{cccc}
\hline
& \multicolumn{3}{c}{Sampling Strategy Costs (\$)} \\
\cline{2-4}
\# of Scenarios    & MC & IS & BQ\\
\hline
5      & \num{1.355e7}   & \num{1.363e7} & \num{1.356e7}      \\
10       & \num{1.358e7}   & \num{1.365e7} & \num{1.356e7}      \\
20       & \num{1.360e7}   & \num{1.366e7} & \num{1.360e7}      \\
50      & \num{1.365e7}   & \num{1.366e7} & \num{1.360e7}      \\
\hline
\end{tabular}
\end{table}

\begin{table}[!t]
\caption{Economic Dispatch Second Stage Costs}
\label{tab:costs3}
\centering
\begin{tabular}{cccc}
\hline
& \multicolumn{3}{c}{Sampling Strategy Costs (\$)} \\
\cline{2-4}
\# of Scenarios    & MC & IS & BQ\\
\hline
5      & \num{9.839e5}   & \num{2.010e5} & \num{3.838e5}      \\
10       & \num{5.101e5}   & \num{1.688e5} & \num{3.838e5}      \\
20       & \num{3.073e5}   & \num{1.732e5} & \num{1.858e5}      \\
50      & \num{2.099e5}   & \num{1.648e5} & \num{1.834e5}      \\
\hline
\end{tabular}
\end{table}

\begin{figure*}[!t]
\centering
\includegraphics[width=0.9\textwidth]{samplingtimeseries.png}
\caption{Costs of the economic dispatch decisions from the three sampling strategies for 5, 10, 20, and 50 scenarios.  While the Bayesian quadrature approach triggers expensive second stage costs more often than importance sampling, it achieves lower total costs for 20 and 50 scenarios with much lower computational costs.}
\label{fig:scenarioTimeseries}
\end{figure*}

With only 5 or 10 scenarios to consider, the importance sampling approach yields the lowest total costs.  We attribute this performance to the additional information available from the deterministic solution.  Recall that an estimate of the loss function $L\left(\bm{x},\bm{\xi}\right)$ is required to determine the importance distribution.  We approximate the loss function by solving a deterministic dispatch problem by assuming a perfect wind forecast.  We then use the trapezoidal rule to calculate the normalizing factor that makes the importance distribution a proper pdf.  Although this is a flawed estimate of the loss function, it is still valuable when only a few scenarios are available for optimization.

For 20 or more scenarios, the Bayesian quadrature approach yields a lower total cost.  At this level of sampling, the Gaussian process approximation of the loss function is more accurate than the deterministic approximation made in importance sampling.  We also note that the Bayesian quadrature approach is much less computationally expensive than importance sampling because it avoids solving an additional deterministic problem and integrating to normalize the importance distribution.  Interestingly, the Bayesian quadrature solution yields a lower total cost than importance sampling despite triggering loss-of-load penalties more often.  This is a reminder that the economically optimal dispatch is not necessarily the most reliable.

\section{Conclusion}
This paper introduced two new techniques, importance sampling and Bayesian quadrature, to create scenarios of wind forecast error events for use in solving stochastic economic dispatch via a two-stage recourse problem.  The proposed techniques were compared against a Monte Carlo sampling strategy in a weeklong study of 5 minute economic dispatch using realistic generation, load, and wind forecast data from the RTS-GMLC test case and NREL's WIND Toolkit dataset.  Both importance sampling and Bayesian quadrature achieve a reduction in estimator variance as compared to the baseline Monte Carlo sampling, which ultimately yield lower total costs for the two-stage economic dispatch decisions.  Furthermore, the scenarios selected by the improved strategies are drawn further out into the tails of the distribution of possible wind ramp events than with Monte Carlo sampling.  The resulting dispatch decisions improve grid reliability and robustness by finding feasible solutions that reduce loss of load, while still maintaining economically optimal dispatch decisions.

We also found that the computational cost of determining sample locations is lower for Bayesian quadrature.  Importance sampling suffers from the \textit{curse of circularity} in that the normalization factor to make the importance distribution a proper probability density function requires knowing the expectation one is trying to estimate in the first place.  To overcome this, we solve a deterministic economic dispatch problem to estimate the loss function and then numerically integrate it to approximate the normalization constant.  This additional solve increases the computation cost of importance sampling relative to Bayesian quadrature.  Furthermore, in the Bayesian quadrature approach the sample points can be computed once and reused at successive timesteps because it is a true quadrature approach rather than a stochastic sampling approach.  The Bayesian quadrature samples only depend on the nominal wind pdf and the Gaussian process kernel, allowing them to be computed offline for wind distributions that are specific to a time of day, season, or current power output.  Our analysis simply used a Student's t-distribution to represent the distribution of wind ramps over the course of a year, however more condition-specific distributions would likely improve performance even further.

The most economic dispatch decisions are found when using importance sampling with ten or fewer scenarios, and with Bayesian quadrature for more than ten scenarios.  We attribute this to the extra information about the loss function that is gained through its deterministic approximation while computing the importance distribution.  With more than ten scenarios, the Bayesian quadrature approach can develop a better loss function surrogate and find a lower cost dispatch strategy.  While both proposed strategies improve the reliability from the baseline, the most economic strategy is not necessarily the most reliable.  The solutions found via Bayesian quadrature for 20 and 50 scenarios produce the lowest total cost as shown in Table I, but as shown in Table III, they involve more recourse or second stage costs than the less economic strategies found via importance sampling, reflecting increased number of expected loss of load events.

Applying these strategies to a realistic test case of 5 minute economic dispatch with time varying loads, stochastic wind generation, and realistic costs yields improvements in total cost and in reliability represented by the second stage recourse costs.  Depending on the number of scenarios used, total costs are reduced by $4.8\%$ - $0.6\%$ relative to the Monte Carlo baseline.  A large component of the improvements in overall economic optimality is actually due to improvements in second stage costs that represent consumption of reserves or loss of load.  The economically optimal dispatch strategy reduces these costs by $79.6\%$ - $12.6\%$, demonstrating significant benefits to grid reliability and reserve requirements.  Employing these advanced scenario generation strategies can significantly reduce costs and improve the reliability of operating the grid with high levels of stochastic renewable energy production.

\section*{Acknowledgment}
This work was supported by the U.S. Department of Energy under Contract No. DE-AC36-08GO28308 with Alliance for Sustainable Energy, LLC, the Manager and Operator of the National Renewable Energy Laboratory. Funding provided by the Department of Energy's Grid Modernization Lab Consortium and Exascale Computing Project.

\bibliography{bibliography}
\bibliographystyle{IEEEtran}

\end{document}

config: arxiv


TITLE: Find the values of $k$ that satisfy $x_1> 1$ and $x_2 <1.$
QUESTION [2 upvotes]: Here I found a question that caught my attention. But then I did not see it again. The question was:
The quadratic equation is given as follows:

$$(2k+1)x^2-kx+k-2=0$$

Find the values of $k$ that satisfy $x_1> 1$ and $x_2 <1.$
My way:
$$x_1=\frac{k+\sqrt{k^2-4(2k+1)(k-2)}}{2(2k+1)} \\
x_2=\frac{k-\sqrt{k^2-4(2k+1)(k-2)}}{2(2k+1)}\\ \begin{cases} x_1>1 & \\ x_2<1 &  \end{cases} \Rightarrow \begin{cases} \frac{k+\sqrt{k^2-4(2k+1)(k-2)}-2(2k+1)}{2(2k+1)} >0 & \\ \frac{k-\sqrt{k^2-4(2k+1)(k-2)}-2(2k+1)}{2(2k+1)}<0  &  \end{cases} \Rightarrow \begin{cases} 2(2k+1)\left(\sqrt{-7k^2+12k+8}-(3k+2) \right)>0 & \\2(2k+1)\left(\sqrt{-7k^2+12k+8}+(3k+2) \right)>0  &  \end{cases} \Rightarrow 4(2k+1)^2×\left( -7k^2+12k+8-(3k+2)^2\right)>0 \Rightarrow -16k^2+4>0 \Rightarrow (2k-1)(2k+1)<0 \Rightarrow k\in \left(-\frac12 ; \frac 12 \right)$$
Actually, I doubt my solution. Can you confirm that the solution is right or wrong?
Thank you!

REPLY [2 votes]: As I've already commented, it seems that you showed that if $x_1\gt 1$ and $x_2\lt 1$, then $k\in\left(-\frac 12,\frac 12\right)$, and that you have not shown that if $k\in\left(-\frac 12,\frac 12\right)$, then $x_1\gt 1$ and $x_2\lt 1$.
This answer shows in your method that $x_1\gt 1$ and $x_2\lt 1$ is equivalent to $k\in\left(-\frac 12,\frac 12\right)$.
$$\small\begin{align}&x_1\gt 1\quad \text{and}\quad x_2\lt 1\\\\&\iff 2k+1\not=0\quad\text{and}\quad (-k)^2-4(2k+1)(k-2)\gt 0\\\\&\qquad\qquad\text{and}\quad \frac{k-\sqrt{(-k)^2-4(2k+1)(k-2)}}{2(2k+1)}\lt 1\quad\text{and}\quad \frac{k+\sqrt{(-k)^2-4(2k+1)(k-2)}}{2(2k+1)}\gt 1
\\\\&\iff k\not=-\frac 12\quad\text{and}\quad \frac{6-2\sqrt{23}}{7}\lt k\lt \frac{6+2\sqrt{23}}{7}\\\\&\qquad\qquad \text{and}\quad \frac{-3k-2-\sqrt{-7k^2+12k+8}}{2k+1}\lt 0\quad\text{and}\quad \frac{-3k-2+\sqrt{-7k^2+12k+8}}{2k+1}\gt 0
\\\\&\iff k\in\left(\frac{6-2\sqrt{23}}{7},-\frac 12\right)\cup\left(-\frac 12,\frac{6+2\sqrt{23}}{7}\right)\\\\&\qquad\qquad \text{and}\quad \frac{-3k-2-\sqrt{-7k^2+12k+8}}{2k+1}\lt 0\quad\text{and}\quad \frac{-3k-2+\sqrt{-7k^2+12k+8}}{2k+1}\gt 0\\\\&\iff \left[k\in\left(\frac{6-2\sqrt{23}}{7},-\frac 12\right)\ \ \text{and}\ \ \sqrt{-7k^2+12k+8}\lt -3k-2\ \ \text{and}\ \ \sqrt{-7k^2+12k+8}\lt 3k+2\right]\\\\&\qquad\qquad\ \text{or}\quad \left[k\in \left(-\frac 12,\frac{6+2\sqrt{23}}{7}\right)\ \ \text{and}\ \ \sqrt{-7k^2+12k+8}\gt -3k-2\ \ \text{and}\ \ \sqrt{-7k^2+12k+8}\gt 3k+2\right]\\\\&\iff k\in \left(-\frac 12,\frac{6+2\sqrt{23}}{7}\right)\quad\text{and}\quad \sqrt{-7k^2+12k+8}\gt |3k+2|\\\\&\iff k\in \left(-\frac 12,\frac{6+2\sqrt{23}}{7}\right)\quad\text{and}\quad -7k^2+12k+8\gt (3k+2)^2\\\\&\iff k\in \left(-\frac 12,\frac{6+2\sqrt{23}}{7}\right)\quad\text{and}\quad k\in\left(-\frac 12,\frac 12\right)\\\\&\iff k\in\left(-\frac 12,\frac 12\right)\qquad\blacksquare\end{align}$$