\section{Introduction}
Diophantine geometry studies rational solutions of the  polynomial equations in terms of
geometry of  the corresponding solutions in complex numbers [Hindry \& Silverman 2000]  \cite[Introduction]{HS}.
The idea that arithmetic depends on geometry is amazingly fruitful,
e.g. the Grothendieck's concept of a scheme over the Dedekind domain  led to the proof of Weil's
Conjectures.    Usually  the arithmetic schemes corresponding to the polynomial rings with rational  coefficients support
morphisms defined over the algebraic closure of the underlying field.
As a result, such morphisms do not  preserve   rational points  on the algebraic varieties.

It is known since the works of J.~-P.~Serre and P.~Gabriel,  that algebraic geometry can be recast
in terms of non-commutative rings $R$ \cite[Chapter 13]{N}.  Unlike commutative rings, the  $R$ are
 usually simple rings.  In particular,  the spectrum of the ring  $R$ is a singleton
 and  one cannot localize $R$ or define a scheme.
However  the simplicity of  $R$ is beneficial  for the diophantine geometry
 and the aim of our note is to show that the  morphisms of $R$ preserve rational points on  the  algebraic variety.
Such a property allows to recover  subtle arithmetic invariants  in terms of the invariants of the ring $R$,
see Section 4.  Let us recall some definitions.

Denote by $V$  a complex projective variety and let $\mathscr{A}_V$ be the
Serre $C^*$-algebra, i.e.  the norm closure of a self-adjoint representation of the twisted
 homogeneous coordinate ring of  $V$ by the bounded linear operators on a Hilbert space
 \cite[Section 5.3.1]{N}.
 Recall that the $C^*$-algebras  $\mathscr{A}$ and $\mathscr{A}'$ are said to be (strongly)  Morita equivalent,
 if $\mathscr{A}\otimes\mathbf{K}\cong \mathscr{A}'\otimes\mathbf{K}$,
 where $\mathbf{K}$ is the $C^*$-algebra of all compact operators on a Hilbert space and
 $\cong$ is the $C^*$-algebra isomorphism [Blackadar 1986]  \cite[13.7.1(c)]{B}.
The Morita equivalence can be viewed as a generalization
 of an isomorphism meaning
 that the $C^*$-algebras $\mathscr{A}$ and $\mathscr{A}'$ have an isomorphic  category
 of the finitely generated projective modules over $\mathscr{A}$ and $\mathscr{A}'$, respectively.
 It is easy to see, that if $\mathscr{A}\cong\mathscr{A}'$ are isomorphic $C^*$-algebras, then they
 are Morita equivalent to each other but not {\it vice versa}.
 It is known, that the Morita equivalences of  the $C^*$-algebra $\mathscr{A}_V$
 preserve the complex points on $V$,  i.e. the projective variety $V$ is $\mathbf{C}$-isomorphic
 to a projective variety $V'$ if and only if the corresponding Serre $C^*$-algebras are Morita
 equivalent \cite[Theorem 5.3.3]{N}.

 Let $k$ be a number field and let $V(k)$ be a projective variety over the field $k$.
 Consider a   $C^*$-algebra $M_n(\mathscr{A}_{V(k)})$
consisting of the $n\times n$ matrices with the entries in $\mathscr{A}_{V(k)}$,
see  [Blackadar 1986]  \cite[p. 16]{B} for the details.
 It follows from the K-theory of $\mathscr{A}_{V(k)}$,  that the $C^*$-algebra $M_n(\mathscr{A}_{V(k)})$ is simple,
 see  lemma \ref{lm3.1}.
 Denote by $K$ an extension of the field $k$ obtained by adjoining
 the roots of all  polynomials of degree $n$ over $k$.
 Our main results can be stated as follows.
\begin{theorem}\label{thm1.1}
 Projective variety $V(k)$ is $K$-isomorphic
 to a projective variety $V'(k)$  if and only if

 \begin{equation}\label{eq1.1}
 M_n\left(\mathscr{A}_{V(k)}\right)\cong  M_n\left(\mathscr{A}_{V'(k)}\right).
 \end{equation}

\end{theorem}
\begin{remark}\label{rmk1.2}
Roughly speaking, theorem \ref{thm1.1} says that the isomorphisms
of the $C^*$-algebra  $M_n(\mathscr{A}_{V(k)})$ preserve the $K$-rational
points on the variety $V(k)$.
\end{remark}
\begin{corollary}\label{cor1.3}
Projective variety $V(k)$ is  $k$-isomorphic ($\mathbf{C}$-isomorphic, resp.) to a projective variety  $V'(k)$
if and only if the Serre $C^*$-algebra  $\mathscr{A}_{V(k)}$ is  isomorphic (Morita equivalent, resp.) to the Serre $C^*$-algebra
  $\mathscr{A}_{V'(k)}$.
\end{corollary}
\begin{remark}\label{rmk1.4}
An independent proof of corollary \ref{cor1.3} using the Galois cohomology
can be found in \cite{Nik1}. Our current proof is more elegant and follows from
(\ref{eq1.1}) and an analog of the Zariski's Lemma for the algebra $\mathscr{A}_{V}$ (lemma \ref{lm3.2}).
\end{remark}
The paper is organized as follows. The preliminary facts can be found in Section 2.
Theorem \ref{thm1.1} and corollary \ref{cor1.3} are proved in Section 3.
An application of  corollary \ref{cor1.3} to the Faltings Theorem is considered in Section 4.

\section{Preliminaries}
We briefly review noncommutative algebraic geometry, operator algebras
and Serre $C^*$-algebras. For a detailed exposition we refer the reader
to [Stafford \& van ~den ~Bergh 2001]  \cite{StaVdb1}, [Blackadar 1986]  \cite[Chapter II]{B}
and \cite[Section 5.3.1]{N}, respectively.

\subsection{Noncommutative algebraic geometry}
Let $V$ be a complex projective variety.
For an automorphism $\sigma: V\to V$ and an invertible
sheaf $\mathcal{L}$ of the linear forms on $V$,  one can construct
a twisted homogeneous coordinate ring $B(V, \mathcal{L}, \sigma)$ of the variety $V$,
i.e. a non-commutative ring such that:
\begin{equation}\label{eq2.1}
Mod~(B(V, \mathcal{L}, \sigma)) ~/~Tors \cong ~Coh~(V),
\end{equation}
where $Mod$ is the category of graded left modules over the graded ring
$B(V, \mathcal{L}, \sigma)$, $Tors$ the full subcategory of $Mod$ of the
torsion modules and $Coh$ the category of quasi-coherent sheaves on the
variety  $V$ [Stafford \& van ~den ~Bergh 2001]  \cite[p.180]{StaVdb1}.
The correspondence $V\mapsto B(V, \mathcal{L}, \sigma)$ is a functor
which maps $\mathbf{C}$-isomorphic projective varieties to the Morita equivalent
rings $B(V, \mathcal{L}, \sigma)$.
\begin{example}\label{exm2.1}
Let $k$ be a field and $U_{\infty}(k)$ the algebra of polynomials
over $k$ generated by  two non-commuting variables $x_1$ and $x_2$ satisfying
the  relation $x_1x_2-x_2x_1-x_1^2=0$.   Let $\mathbb{P}^1(k)$ be the projective
line over $k$.  Then $B(V, \mathcal{L}, \sigma)\cong U_{\infty}(k)$ and $V\cong \mathbb{P}^1(k)$
satisfy  (\ref{eq2.1}).   Notice  that the $B(V, \mathcal{L}, \sigma)$ is far from being a commutative ring.
\end{example}
\begin{example}\label{exm2.2}
Let $S(\alpha,\beta,\gamma)$ be the  Sklyanin algebra, i.e
a free $\mathbf{C}$-algebra on  four generators
$\{x_1, \dots, x_4\}$ satisfying six quadratic relations:
$x_1x_2-x_2x_1 = \alpha(x_3x_4+x_4x_3)$,
\quad $x_1x_2+x_2x_1 = x_3x_4-x_4x_3$,
~$x_1x_3-x_3x_1 = \beta(x_4x_2+x_2x_4)$,
~$x_1x_3+x_3x_1 = x_4x_2-x_2x_4$,
~$x_1x_4-x_4x_1 = \gamma(x_2x_3+x_3x_2)$ and
$x_1x_4+x_4x_1 = x_2x_3-x_3x_2$,
where $\alpha,\beta,\gamma\in \mathbf{C}$ and
$\alpha+\beta+\gamma+\alpha\beta\gamma=0$.
Let  $\mathcal{E}\subset \mathbf{C}P^3$
be an elliptic curve given in the Jacobi form,  i.e.  as  an  intersection of two quartics:
\begin{equation}\label{eq2.2}
\left\{
\begin{array}{ccc}
u^2+v^2+w^2+z^2 &=&  0,\\
{1-\alpha\over 1+\beta}v^2+{1+\alpha\over 1-\gamma}w^2+z^2  &=&  0.
\end{array}
\right.
\end{equation}
 Then $B(V, \mathcal{L}, \sigma)\cong S(\alpha,\beta,\gamma)$ and $V\cong \mathcal{E}$
satisfy  (\ref{eq2.1}). In other words,   the   $S(\alpha,\beta,\gamma)$ is a twisted homogenous
coordinate ring of the elliptic curve (\ref{eq2.2}).
\end{example}

\subsection{Operator algebras}
The $C^*$-algebra is an algebra  $\mathscr{A}$ over $\mathbf{C}$ with a norm
$a\mapsto ||a||$ and an involution $\{a\mapsto a^* ~|~ a\in \mathscr{A}\}$  such that $\mathscr{A}$ is
complete with  respect to the norm, and such that $||ab||\le ||a||~||b||$ and $||a^*a||=||a||^2$ for every  $a,b\in \mathscr{A}$.
Each commutative $C^*$-algebra is  isomorphic
to the algebra $C_0(X)$ of continuous complex-valued
functions on some locally compact Hausdorff space $X$.
Any other  algebra $\mathscr{A}$ can be thought of as  a noncommutative
topological space.
\begin{example}\label{exm2.3}
An $n$-dimensional   noncommutative  torus  is the universal  $C^*$-algebra $\mathscr{A}_{\Theta_n}$
generated by $n$ unitary operators $U_1,\dots, U_n$ satisfying the commutation relations
$\{U_jU_i=\exp ~(2\pi i ~\theta_{ij}) ~U_iU_j ~|~1 \le i,j\le n\}$ which depend on a skew-symmetric matrix
\begin{equation}\label{skew}
\Theta_n=\left(
\begin{matrix}
              0 & \theta_{12}  & \dots & \theta_{1n}\cr
             -\theta_{12} & 0  & \dots & \theta_{2n}\cr
              \vdots         & \vdots         & \ddots   &\vdots\cr
             -\theta_{1n} & -\theta_{2n} & \dots & 0
             \end{matrix}
              \right),
              \quad\hbox{where}\quad  \theta_{ij}\in \mathbf{R}.
              \end{equation}
\end{example}

\bigskip
By $M_{\infty}(\mathscr{A})$
one understands the algebraic direct limit of the $C^*$-algebras
$M_n(\mathscr{A})$ under the embeddings $a\mapsto ~\mathbf{diag} (a,0)$.
The direct limit $M_{\infty}(\mathscr{A})$  can be thought of as the $C^*$-algebra
of infinite-dimensional matrices whose entries are all zero except for a finite number of the
non-zero entries taken from the $C^*$-algebra $\mathscr{A}$.
Two projections $p,q\in M_{\infty}(\mathscr{A})$ are equivalent, if there exists
an element $v\in M_{\infty}(\mathscr{A})$,  such that $p=v^*v$ and $q=vv^*$.
The equivalence class of projection $p$ is denoted by $[p]$.
We write $V(\mathscr{A})$ to denote all equivalence classes of
projections in the $C^*$-algebra $M_{\infty}(\mathscr{A})$, i.e.
$V(\mathscr{A}):=\{[p] ~:~ p=p^*=p^2\in M_{\infty}(\mathscr{A})\}$.
The set $V(\mathscr{A})$ has the natural structure of an abelian
semi-group with the addition operation defined by the formula
$[p]+[q]:=\mathbf{diag}(p,q)=[p'\oplus q']$, where $p'\sim p, ~q'\sim q$
and $p'\perp q'$.  The identity of the semi-group $V(\mathscr{A})$
is given by $[0]$, where $0$ is the zero projection.
By the $K_0$-group $K_0(\mathscr{A})$ of the unital $C^*$-algebra $\mathscr{A}$
one understands the Grothendieck group of the abelian semi-group
$V(\mathscr{A})$, i.e. a completion of $V(\mathscr{A})$ by the formal elements
$[p]-[q]$.
\begin{example}\label{exm2.4}
{\bf ([Rieffel 1990] \cite{Rie1})}
\begin{equation}
K_0(\mathscr{A}_{\Theta_n})\cong \mathbf{Z}^{2^{n-1}}.
\end{equation}
\end{example}

\subsection{Serre $C^*$-algebra}
Let $V$ be a complex projective variety and let  $B(V, \mathcal{L}, \sigma)$ be
 twisted homogeneous coordinate ring of $V$ (Section 2.1).  Denote by $\mathscr{H}$
 a Hilbert space and let $\mathbf{B}(\mathscr{H})$ be the algebra of bounded linear operators
 on $\mathscr{H}$.  Fix a self-adjoint (i.e. $\ast$-invariant) representation

 \begin{equation}\label{eq2.4}
 \rho: B(V, \mathcal{L}, \sigma)\longrightarrow \mathbf{B}(\mathscr{H}).
 \end{equation}

By the {\it Serre $C^*$-algebra}  of $V$ one understands the closure
$\mathscr{A}_V$ of the $\ast$-algebra $\rho(B(V, \mathcal{L}, \sigma))$
in the norm topology on the space $\mathbf{B}(\mathscr{H})$.
\begin{example}\label{exm2.5}
Let $V\cong \mathbf{A}_g$ be a $g$-dimensional abelian variety,
where $g\ge 1$.  The corresponding Serre $C^*$-algebra
$\mathscr{A}_V\cong \mathscr{A}_{\Theta_{2g}}$ is a $2g$-dimensional
noncommutative torus, see  Example \ref{exm2.3}.
In particular for $g=1$ and $V\cong\mathcal{E}$ being an elliptic curve
as in Example \ref{exm2.2},  one gets  $\mathscr{A}_V\cong \mathscr{A}_{\theta}$
is an  irrational rotation algebra, i.e. the universal $C^*$-algebra on two generators $U_1$ and $U_2$
satisfying the relation $U_2U_1=e^{2\pi i\theta}U_1U_2$ for a constant $\theta\in\mathbf{R}$.
This can be  proved by comparing the generators
and relations of  the Sklyanin algebra $S(\alpha,\beta,\gamma)$ given in Example \ref{exm2.2}
with such for the algebra   $\mathscr{A}_{\theta}$ \cite[Section 1.3]{N}.
\end{example}

\section{Proofs}
\subsection{Proof of theorem \ref{thm1.1}}
We split the proof in a series of  lemmas.

\begin{lemma}\label{lm3.1}
The  $M_n\left(\mathscr{A}_{V(k)}\right)$ is
a simple $C^*$-algebra for each $n\ge 1$.
\end{lemma}
\begin{proof}
(i) Let us prove  that the $\mathscr{A}_{V(k)}$ is a simple $C^*$-algebra.
Indeed,  any Serre $C^*$-algebra $\mathscr{A}_{V}$ is a crossed product $C^*$-algebra \cite[Theorem 5.3.4]{N}.
Denote by $\mathbb{A}$ an Approximately Finite-dimensional (AF-) $C^*$-algebra \cite[Section 3.5]{N},
such that   $\mathscr{A}_{V}\hookrightarrow \mathbb{A}$ is an embedding  and  $K_0(\mathscr{A}_{V})\cong K_0 (\mathbb{A})$.
In particular, the $\mathscr{A}_{V(k)}$ embeds into an AF-algebra $\mathbb{A}$  of stationary type \cite[Section 3.5.2]{N}.
Such AF-algebras are always simple,  since the corresponding dimension group $(K_0 (\mathbb{A}), K_0^+ (\mathbb{A}))$
\cite[Definition 3.5.2]{N} does not have order-ideals.  Since $K_0^+(\mathscr{A}_{V(k)})\cong K_0^+ (\mathbb{A})$,
we conclude that the Serre $C^*$-algebra  $\mathscr{A}_{V(k)}$ is  simple.

\medskip
(ii)  Let us prove  that the  $M_n\left(\mathscr{A}_{V(k)}\right)$ is a simple $C^*$-algebra.
Recall that there exists a one-to-one correspondence between the two-sided ideals of the
matrix algebra $M_n\left(\mathscr{A}_{V(k)}\right)$ and the two-sided of  $\mathscr{A}_{V(k)}$.
In particular, the  $M_n\left(\mathscr{A}_{V(k)}\right)$ is simple if and only if the  $\mathscr{A}_{V(k)}$
is simple. It follows from item (i), that the matrix $C^*$-algebra $M_n\left(\mathscr{A}_{V(k)}\right)$
is simple for each $n\ge 1$.  Lemma \ref{lm3.1} is proved.
\end{proof}

\begin{lemma}\label{lm3.2}
{\bf (Zariski's Lemma for the algebra $\mathscr{A}_V$)}
 If  $K$ is an extension of the number field $k$ obtained  by adjoining
 the roots of all  polynomials of degree $n$ over $k$,  then

 \begin{equation}
 \mathscr{A}_{V(K)}\cong M_n\left(\mathscr{A}_{V(k)}\right).
 \end{equation}

\end{lemma}
\begin{proof}
(i)
Recall that the $\mathscr{A}_{V(k)}$  is the norm-closure of
a self-adjoint representation of   the $k$-algebra
$B:=B(V(k), \mathcal{L}, \sigma)$ of  non-commutative polynomials
 over the number field $k$,  see   Example \ref{exm2.2} and formula (\ref{eq2.4}).  Consider
 the Frobenius companion matrix given by an element
\begin{equation}\label{eq3.2}
\mathscr{C}=
\left(
\begin{matrix}
0 & 0 & \dots & 0 &-c_0\cr
1 & 0 & \dots &  0 &-c_1\cr
\vdots &\vdots & & \vdots& \vdots\cr
0 & 0 & \dots & 1 & -c_{n-1}
\end{matrix}
\right),
\qquad c_i\in k
\end{equation}
of the $K$-algebra $M_n (B)$.
To determine the field of constants $K$ of the algebra  $M_n (B)$,
denote by  $I$ the identity of  $M_n (B)$
and consider the set $\{c I ~|~c\in K\}$.  It is easy to see, that
\begin{equation}\label{eq3.3}
c I=\mathscr{C}\in   M_n (B)
\end{equation}
 if an only if
$c$ is a root of the characteristic polynomial of matrix $\mathscr{C}$
given by (\ref{eq3.2}), i.e.
\begin{equation}\label{eq3.4}
Char ~\mathscr{C}=x^n+c_{n-1}x^{n-1}+\dots+c_0.
\end{equation}

We conclude from (\ref{eq3.3})  that the field of constants $K$ of the algebra $M_n (B)$
must contain all roots of the  polynomial (\ref{eq3.4}). Since $c_i\in k$, these  roots are algebraic numbers of degree $n$ over $k$.

One can run through the set of all matrices $\mathscr{C}$ given by formula (\ref{eq3.2}) with $c_0\ne 0$.
It is clear that in this way we get all algebraic extensions of degree $n$ over the number field $k$.
Thus the field of constants $K$ of the algebra   $M_n (B)$
is obtained from $k$ by adjoining  all algebraic numbers of degree $n$ over $k$.

\bigskip
(ii) It follows from item (i) that the  $M_n (B)$ is an algebra of non-commutative polynomials
over the field $K$.  Since the field of constants of the twisted homogeneous coordinate ring
coincides with the field of definition of projective variety (see Example \ref{exm2.2}),
we conclude that

\begin{equation}\label{eq3.5}
M_n(B)\cong B(V(K), ~\mathcal{L}, ~\sigma).
\end{equation}

\bigskip
(iii) It remains to compare (\ref{eq3.5}) with the notation $B\cong B(V(k), ~\mathcal{L}, ~\sigma)$
and  to obtain a ring  isomorphism:

\begin{equation}\label{eq3.6}
B(V(K), ~\mathcal{L}, ~\sigma)\cong
M_n(B(V(k), ~\mathcal{L}, ~\sigma)).
\end{equation}
Taking the norm-closure of a self-adjoint representation (\ref{eq2.4})
 of the rings  in the both sides of (\ref{eq3.6}),  one gets a $C^*$-algebra
 isomorphism:

\begin{equation}\label{eq3.7}
\mathscr{A}_{V(K)}\cong M_n\left(\mathscr{A}_{V(k)}\right).
\end{equation}

Lemma \ref{lm3.2} is proved.
\end{proof}

\begin{lemma}\label{lm3.3}
 The varieties $V(k)$ and $V'(k)$ are $K$-isomorphic
if and only if
\begin{equation}
 M_n\left(\mathscr{A}_{V(k)}\right)\cong  M_n\left(\mathscr{A}_{V'(k)}\right).
\end{equation}
\end{lemma}
\begin{proof}
(i) Let $V(k)\cong V'(k)$ be a pair of isomorphic projective varieties defined over the
number field $k$. The canonical embedding
\begin{equation}
\mathscr{A}_{V(k)}\hookrightarrow M_n\left(\mathscr{A}_{V(k)}\right)\cong \mathscr{A}_{V(K)}
\end{equation}
gives rise to an embedding $i$ of the variety $V(k)$ into the variety $V(K)$.
Consider a commutative diagram in Figure 1, where $\psi$ is an isomorphism defined over
the number field $K$.
Since $\varphi=i\circ\psi\circ  i^{-1}$ is a composition of  maps  defined over the number field  $K$,  we conclude that the isomorphism
$\varphi: V(k)\to V'(k)$ must also be  defined over the  field $K$.
\begin{figure
\begin{picture}(300,110)(-60,-5)
\put(30,70){\vector(0,-1){35}}
\put(140,70){\vector(0,-1){35}}
\put(55,23){\vector(1,0){60}}
\put(55,83){\vector(1,0){60}}
\put(15,20){$V'(k)$}
\put(128,20){$V'(K)$}
\put(17,80){$V(k)$}
\put(125,80){$V(K)$}
\put(80,30){$i$}
\put(80,90){$i$}
\put(10,50){$\varphi$}
\put(125,50){$\psi$}
\end{picture}
\caption{}
\end{figure}

\bigskip
(ii) Let $\psi$ be the $K$-isomorphism between the varieties $V(K)$ and $V'(K)$ shown on the diagram in
Figure 1.  The $\psi $ gives rise to an isomorphism  $\mathscr{A}_{V(K)}\cong \mathscr{A}_{V'(K)}$
of the corresponding Serre $C^*$-algebras.

In view of lemma \ref{lm3.2},  one gets the isomorphisms $\mathscr{A}_{V(K)}\cong M_n\left(\mathscr{A}_{V(k)}\right)$ and
$\mathscr{A}_{V'(K)}\cong M_n\left(\mathscr{A}_{V'(k)}\right)$.
Comparing the latter with the results of item (i), we conclude that
the varieties $V(k)\cong V'(k)$ are $K$-isomorphic
if and only if
\begin{equation}
 M_n\left(\mathscr{A}_{V(k)}\right)\cong  M_n\left(\mathscr{A}_{V'(k)}\right).
\end{equation}

Lemma \ref{lm3.3} is proved.
\end{proof}

\bigskip
Theorem \ref{thm1.1} follows from lemma \ref{lm3.3}.

\subsection{Proof of corollary \ref{cor1.3}}
(i) Let us prove that
the projective variety $V(k)$ is  $k$-isomorphic  to a projective variety  $V'(k)$
if and only if the Serre $C^*$-algebra  $\mathscr{A}_{V(k)}$ is  isomorphic  to the Serre $C^*$-algebra
  $\mathscr{A}_{V'(k)}$.  Indeed, suppose that in theorem \ref{thm1.1} we have  $n=1$.
  In this case the field $K$ is an algebraic extension of the field $k$ of degree $1$.
  In other words, the fields $k\cong K$ are isomorphic.

  On the other hand, theorem \ref{thm1.1} says that the variety $V(k)$ is $k$-isomorphic  to the variety
  $V'(k)$ if and only if    $\mathscr{A}_{V(k)}\cong \mathscr{A}_{V'(k)}$ by formula
  (\ref{eq1.1}) with $n=1$.  Item (i) is proved.

\bigskip
(ii) Let us prove that
the projective variety $V(k)$ is  $\mathbf{C}$-isomorphic to a projective variety  $V'(k)$
if and only if the Serre $C^*$-algebra  $\mathscr{A}_{V(k)}$ is  Morita equivalent  to the Serre $C^*$-algebra
  $\mathscr{A}_{V'(k)}$.

We shall denote by $\bar k$ the algebraic closure of the number field $k$ obtained by adjoining the roots of
all polynomials of finite degree over $k$.  It is easy to see, that  this case of theorem \ref{thm1.1} corresponds to $n=\infty$ and
$K\cong \bar k$.

Recall that the $C^*$-algebra $M_n\left(\mathscr{A}_{V(k)}\right)$ can be written in the form of a tensor product, i.e.
\begin{equation}\label{eq3.11}
M_n\left(\mathscr{A}_{V(k)}\right)\cong \mathscr{A}_{V(k)}\otimes M_n(\mathbf{C}).
\end{equation}

 On the other hand, the canonical inclusions of the $C^*$-algebras $\mathbf{C}\subset M_2(\mathbf{C})\subset M_3(\mathbf{C})\subset\dots$
 are known to converge to an inductive limit:
\begin{equation}
\lim_{n\to\infty} M_n(\mathbf{C})\cong \mathbf{K},
\end{equation}
where $\mathbf{K}$ is the $C^*$-algebra of all compact operators on a Hilbert space.
Thus in view of (\ref{eq3.11}), the isomorphism (\ref{eq1.1}) for $n=\infty$ can be written in the form:
\begin{equation}\label{eq3.13}
\mathscr{A}_{V(k)}\otimes\mathbf{K} \cong
\mathscr{A}_{V'(k)}\otimes\mathbf{K}.
\end{equation}
By definition, the isomorphism (\ref{eq3.13}) says that the Serre $C^*$-algebras   $\mathscr{A}_{V(k)}$
and  $\mathscr{A}_{V'(k)}$  are Morita equivalent.

\medskip
It remains to apply the Lefschetz Principle to the number field $\bar k$.  Namely, one can pass from the algebraically
closed field $\bar k$  of characteristic zero to the field of complex numbers $\mathbf{C}$.
This remark finishes the proof of item (ii).

\bigskip
(iii) Corollary \ref{cor1.3} follows from items (i) and (ii).

\section{Faltings Theorem}
In this section we consider an application of the corollary \ref{cor1.3}
to  an alternative proof of the Faltings Theorem about finiteness of the
rational points  on the higher genus curves  [Hindry \& Silverman 2000]  \cite[Theorem E.0.1]{HS}.
Namely, we show that corollary \ref{cor1.3} and a classification of the even-dimensional
noncommutative tori up to an isomorphism (lemma \ref{lm4.2})  imply  the Shafarevich Conjecture for the abelian varieties [Par\v{s}in 1970] \cite[Conjecture S1]{Par1}.
The Faltings Theorem is known to follow from the Shafarevich Conjecture using the  Parshin's Trick.  Let us pass to a detailed argument.

\begin{theorem}\label{thm4.1}
{\bf (Shafarevich Conjecture)}
There exists only a finite number of the abelian varieties $\mathbf{A}_g$ over a number field $k$, such that:

\medskip
(i)   $\mathbf{A}_g$ have a  polarization of  degree $d\ge 1$;

\smallskip
(ii) $\mathbf{A}_g$ have a good reduction outside a fixed finite set $S$ of places of the field $k$.
\end{theorem}
\begin{proof}
We split the proof in two lemmas having an independent interest.

\begin{lemma}\label{lm4.2}
There exists only a finite number of pairwise distinct skew-symmetric matrices $\Theta_{2g}$ given by the formula (\ref{skew}),
such that the corresponding noncommutative tori $\mathscr{A}_{\Theta_{2g}}$ are isomorphic to each other.
\end{lemma}
\begin{proof}
(i) Denote by $\rho_{ij}=\exp~(2\pi i\theta_{ij})$, see Example \ref{exm2.3}.
Then the commutation relations for the noncommutative torus
$\mathscr{A}_{\Theta_{2g}}$ can be written in the form:
\begin{equation}
U_iU_jU_i^{-1}U_j^{-1}=\rho_{ij}\in\mathbf{T},
\end{equation}
where $\mathbf{T}$ is the unit circle.  Since the $n$-dimensional
noncommutative torus depends on ${n(n-1)\over 2}$ real parameters
$\theta_{ij}$, one can identify the $\mathscr{A}_{\Theta_{2g}}$ with the
family of noncommutative tori
\begin{equation}\label{eq4.2}
\left\{\mathscr{A}_{\rho} ~|~ \rho\in\mathbf{T}^{g(2g-1)}\right\}.
\end{equation}

\bigskip
(ii)
 For a dense subalgebra $\mathscr{A}_{\rho}^0\subset \mathscr{A}_{\rho}$ consisting of the polynomials
 in the variables $\{U_k ~|~ 1\le k\le 2g\}$ and an element $x\in \mathscr{A}_{\rho}^0$
 consider a function given by the formula

\begin{equation}\label{eq4.3}
\rho\mapsto ||x||.
\end{equation}
It is easy to see,  that (\ref{eq4.3}) is a (real) analytic function on $\mathbf{T}^{g(2g-1)}$.

\bigskip
(iii)  Lemma \ref{lm4.2} can be proved by contradiction.
Let  $\{\rho_k\}\in \mathbf{T}^{g(2g-1)}$ be an infinite sequence,  such that the $\mathscr{A}_{\rho_k}$ in
 (\ref{eq4.2})  are pairwise isomorphic noncommutative tori.
Since $\mathbf{T}^{g(2g-1)}$ is a compact space, one can always restrict to a Cauchy subsequence
$\{\rho_k'\}$ convergent to a point $\rho\in  \mathbf{T}^{g(2g-1)}$.
Since function (\ref{eq4.3}) is constant on $\rho_k'$ and analytic on an interval $[\rho_k',\rho]$,
we conclude that (\ref{eq4.3}) is constant everywhere on the interval  $[\rho_k',\rho]$.
This is obviously false,  since the interval $[\rho_k',\rho]\subset  \mathbf{T}^{g(2g-1)}$
contains  a continuum of pairwise non-isomorphic noncommutative tori $\mathscr{A}_{\rho}$
 [Rieffel 1990] \cite{Rie1}.  This contradiction finishes the proof of lemma \ref{lm4.2}.
\end{proof}

\bigskip
\begin{example}\label{exm4.3}
Let $g=1$ in lemma \ref{lm4.2}. Then matrix (\ref{skew}) can be written as:
\begin{equation}
\Theta=\left(
\begin{matrix}
0 & \theta\cr
-\theta & 0
\end{matrix}
\right),  \qquad \hbox{where} \quad \theta\in\mathbf{R}.
\end{equation}
It is known, that the noncommutative tori $\mathscr{A}_{\theta}$ and $\mathscr{A}_{\theta'}$
are isomorphic if and only if $\theta'=\pm\theta\mod\mathbf{Z}$ [Rieffel 1990] \cite[p. 200]{Rie1}.
In other words, there are only two distinct values $\rho_1,\rho_2\in\mathbf{T}$, such that
$\rho_{1,2}=\exp (\pm 2\pi i\theta)$ correspond to the isomorphic noncommutative tori.
\end{example}
\begin{remark}
For $g\ge 2$ an explicit formula for the matrices $\Theta_{2g}$ corresponding to the pairwise
isomorphic noncommutative tori is an open problem  [Rieffel 1990] \cite[Question 4.1]{Rie1}.
Lemma \ref{lm4.2} can be viewed as a partial answer to this problem.
\end{remark}

\begin{lemma}\label {lm4.5}
The number of abelian varieties  $\mathbf{A}_g$ over a number field $k$ satisfying conditions
(i) and (ii) of Theorem \ref{thm4.1} is at most finite.
\end{lemma}
\begin{proof}
(i) Let $\mathbf{A}_g$ be a $g$-dimensional abelian variety defined over the number field $k$.
Consider a $k$-isomorphism $\varphi:  \mathbf{A}_g\to \mathbf{A}_g'$ from  $\mathbf{A}_g$
to an abelian variety  $\mathbf{A}_g'$.  The corollary \ref{cor1.3} says that the Serre $C^*$-algebras
$\mathscr{A}_{\mathbf{A}_g}\cong\mathscr{A}_{\Theta_{2g}}$ and  $\mathscr{A}_{\mathbf{A}_g'}\cong\mathscr{A}_{\Theta_{2g}'}$
 must be isomorphic  noncommutative tori, see    Example \ref{exm2.5}.  The corresponding commutative
diagram is shown in Figure 2.

\begin{figure
\begin{picture}(300,110)(-60,-5)
\put(30,70){\vector(0,-1){35}}
\put(140,70){\vector(0,-1){35}}
\put(55,23){\vector(1,0){60}}
\put(55,83){\vector(1,0){60}}
\put(20,20){$\mathscr{A}_{\Theta_{2g}}$}
\put(128,20){$\mathscr{A}_{\Theta_{2g}'}$}
\put(25,80){$\mathbf{A}_g$}
\put(130,80){$\mathbf{A}_g'$}
\put(80,30){$\cong$}
\put(80,90){$\varphi$}
\end{picture}
\caption{}
\end{figure}

\bigskip
(ii)
Let $(\mathbf{A}_g, d, S)$ be the abelian variety $\mathbf{A}_g$ over the field $k$ with a polarization
of degree $d\ge 1$ and good reduction outside a finite set $S$ of places of $k$.
Since $\varphi$ is a $k$-isomorphism of $\mathbf{A}_g$, it must preserve the degree $d$
of polarization and the set $S$ of places of a bad reduction of $\mathbf{A}_g$.
In other words,
\begin{equation}\label{eq4.5}
\varphi(\mathbf{A}_g, d, S)=(\mathbf{A}_g', d, S).
\end{equation}

\bigskip
(iii) Lemma \ref{lm4.2} implies  that there exists at most  finite number of isomorphic
abelian varieties  $\mathbf{A}_g$ satisfying (\ref{eq4.5}).  Indeed, for otherwise
one gets an infinite number of isomorphic noncommutative tori  $\mathscr{A}_{\Theta_{2g}}$,
see Figure 2.  This contradicts lemma \ref{lm4.2} and finishes the proof of
lemma \ref{lm4.5}.
\end{proof}

\bigskip
Theorem \ref{thm4.1} follows from lemma \ref{lm4.5}.
\end{proof}

\begin{corollary}\label{cor4.6}
{\bf (Faltings Theorem)}
If $\mathscr{C}(k)$ is a curve of genus $g\ge 2$, then the set of its $k$-points
is finite.
\end{corollary}
\begin{proof}
The argument based on the Shafarevich Conjecture \ref{thm4.1} is well known.
For the sake of completeness, we repeat it in below.

\bigskip
(i) One can always assume that there exists at least one $k$-point of $\mathscr{C}(k)$.
The Jacobian $Jac~\mathscr{C}$ of $\mathscr{C}(k)$ is an  abelian
variety over the field $k$.  Conversely,  given a principal polarization $d=1$
on  the Jacobian $Jac~\mathscr{C}$, one can recover a curve $\mathscr{C}(k)$
having the same set  $S$ as the abelian variety  $Jac~\mathscr{C}$.
In view of Theorem \ref{thm4.1}, there exists only a finite number of
pairs  $(\mathscr{C}(k), S)$ consisting of the isomorphic curves $\mathscr{C}(k)$
with a fixed set $S$.

\bigskip
(ii) The rest of the proof is a Parshin's Trick described in  [Par\v{s}in 1970] \cite[Th\'eor\`eme 1]{Par1}.
Namely, if $P\in\mathscr{C}(k)$ is a $k$-point,  then there exists a branched covering $\mathscr{C}_P(K)$
of  $\mathscr{C}(k)$, such that the number field $K$ is ramified in $S$ over $k$.
The construction does not depend on the choice of the $k$-point $P$.
It is known from item (i)  that there exists only a finite number of the pairs $(\mathscr{C}_P(K), S)$
and, therefore, only a finite number of the $k$-points on the curve  $\mathscr{C}(k)$.

\bigskip
Corollary \ref{cor4.6} follows from item (ii).
\end{proof}

\bigskip
\begin{remark}
It follows from the proof of theorem \ref{thm4.1} and corollary \ref{cor4.6},
that the number  $|\mathscr{C}(k)|$ of the rational points on the higher genus
curve depends on  the number
of isomorphic $n$-dimensional noncommutative tori  given  by the distinct skew-symmetric
matrices (\ref{skew})  [Rieffel 1990] \cite[Question 4.1]{Rie1}.
To the best of our knowledge, the Rieffel's Problem is unsolved except for  the special cases  (see  example \ref{exm4.3}).
\end{remark}

\bibliographystyle{amsplain}


\section*{Nomenclature}

\subsection{Sets and Indices}
\addcontentsline{toc}{subsection}{Sets and Indices}
\begin{IEEEdescription}[\IEEEusemathlabelsep\IEEEsetlabelwidth{$n\in {\cal{N}}^{\mathrm{T/D}}aa$}]
\item[$\Xi^{[\cdot]}$]{Sets of  variables, where $[\cdot]$ denotes the  regulator (R),   utility (U), wholesale market (W), master (MP), subproblem (SP) and dual wholesale market problems (DW)}
\item[$e\in {\cal{E}}$]{Set of representative operating days}
\item[$i\in {\cal{I}}$]{Set of existing generators}
\item[$i\in \hat{\cal{I}}$]{Set of candidate generators}
\item[${\cal{I}^{\mathrm{R}}}, {\cal{I}^{\mathrm{C}}} \subset {\cal{I}}$]{Set of renewable/controllable generators}
\item[${\cal{I}}_s \subset {\cal{I}}$]{Set of generators in state  $s$}
\item[$l\in {\cal{L}}$]{Set of transmission lines}
\item[$n\in {\cal{N}}$]{Set of transmission nodes}
\item[${\cal{N}}_s \subset {\cal{N}}$]{Set of transmission nodes in state $s$}
\item[$s\in {\cal{S}}$]{Set of states and state regulators}
\item[$t\in {\cal{T}}$]{Set of time intervals}
\item[$ {\cal{T}}^{\mathrm{on/off}}$]{Set of on- and off-peak time intervals}
\item[${r(l), o(l)}$]{Receiving/sending nodes of line $l$}
\item[${n(i)}$]{Node where generator $i$ is located}
\item[${s(i)}$]{State where generator $i$ is located}
\end{IEEEdescription}

\subsection{Parameters}
\addcontentsline{toc}{subsection}{Parameters}
\begin{IEEEdescription}[\IEEEusemathlabelsep\IEEEsetlabelwidth{${G}^{\mathrm{min}/\mathrm{max}}_{it}$\!\!\!\!\!\!\!\!}]
\item[$\eta$]{Tolerance of constraint violations}
\item[$\kappa_{s}$]{Renewable portfolio standard goal of regulator $s$}
\item[$\omega_{e}$]{Probability of operating day $e$}
\item[$\rho_{ite}$]{Forecast factor of generator $i$}
\item[$\sigma_{ite}$]{Normalized standard deviation of the generation forecast error of unit  $i$ [$\mathrm{MW}$]}
\item[$\upsilon_{ite}$]{Normalized mean of the generation forecast error of unit $i$ [$\mathrm{MW}$]}
\item[$\Gamma_i$]{Minimum power output factor of controllable generator $i$, i.e. $0 \leq \Gamma_i\leq 1$ }
\item[$B^{\mathrm{P}}_s$]{Budget for renewable policies of state  $s$}
\item[$C^{\mathrm{inv}}_{i}$]{Capital cost of generator $i$ (prorated on a daily basis using the net present value approach) [$\mathrm{\$/MW}$]}
\item[$C^{\mathrm{g}}_{i}$]{Incremental cost of generator $i$ [$\mathrm{\$/MWh}$]}
\item[${D}_{nte}$]{Real power demand of node $n$ [$\mathrm{MW}$]}
\item[${F}^{\mathrm{max}}_l$]{Apparent flow limit of line $l$ [$\mathrm{MVA}$]}
\item[${G}^{\mathrm{max}\!/\!\mathrm{min}}_{it}$]{Power output limits of existing generator $i$ [$\mathrm{MW}$]}
\item[${H}^{\mathrm{max}\!/\!\mathrm{min}}_i$]{Up/down-ward ramping limits of generator $i$}
\item[${P}_n^{\downarrow,\mathrm{max}}$]{Apparent flow limit of the interface line into node $n$ from the transmission network [$\mathrm{MVA}$]}
\item[$X_l$]{Reactance of transmission line $l$ [$\Omega$]}
\end{IEEEdescription}

\subsection{Variables}
\addcontentsline{toc}{subsection}{Variables}
\begin{IEEEdescription}[\IEEEusemathlabelsep\IEEEsetlabelwidth{$m^{(i)}_s$}]
\item[$\alpha_{it}$]{Participation factor of controllable generator $i$}
\item[$\bm{\epsilon}_{ite}$]{Forecast error of renewable generator $i$}
\item[$\theta_{nte}$]{Voltage angle of transmission node $n$}
\item[$\lambda_{nte}$]{Locational marginal price in node $n$ [$\$/\mathrm{MW}$]}
\item[$\pi_{nt}$]{Retail electricity tariff of node $n$}
\item[$ \pi_n^{\mathrm{on/off}}$]{On- and off-peak retail electricity tariff of node $n$}
\item[$\tau^{\mathrm{c}}$]{Capacity-based renewable energy incentive [$\$/\mathrm{kW}$]}
\item[$\tau^{\mathrm{e}}$]{Energy-based renewable energy incentive [$\$/\mathrm{MWh}$]}
\item[$d_{nte}$]{Flexible real power demand of node $n$ [$\mathrm{MW}$]}
\item[$f_{lte}$]{Real power flow of line $l$ [$\mathrm{MW}$]}
\item[$\bm{g}_{ite}$]{Generation output of generator $i$ [$\mathrm{MW}$]}
\item[$\overline{g}_{ite}$]{Forecast power output (offer into the wholesale market) of generator $i$ [$\mathrm{MW}$]}
\item[${g}_{ite}$]{Dispatch decisions of the wholesale market for generator $i$ [$\mathrm{MW}$]}
\item[$g^{\mathrm{max}}_{i}$]{Capacity of  generator $i$ to be built [$\mathrm{MW}$]}
\item[$p^{\downarrow}_{nte}$]{Interface power flow into node $n$ [$\mathrm{MW}$]}

\end{IEEEdescription}
{ We define random variables in \textbf{bold} and dual variables $\xi_{lts}, \underline{\gamma}_{ite}, \overline{\gamma}_{ite}, \underline{\delta}_{ite}, \overline{\delta}_{ite}$ of constraints in \eqref{LLproblem} in parentheses.}

\section{Introduction}\label{Sec:intro}

The \textit{climate crisis} is the foremost priority of environmental policies around the globe, and reducing $\mathrm{CO_2}$ emissions from the electricity sector is a second-to-none curtailment instrument.
Electricity production accounts for 33\% of annual $\mathrm{CO_2}$ emissions in the U.S. and 65\% of that, i.e. 1,150 MMmt out of 1,763 MMmt, is produced by coal generators \cite{EIA_CO2}.
Thus, many policy-makers aim   to decarbonize electricity production, relying on  a variety of policy instruments, including forced power plant retirements, carbon taxes, cap-and-trade programs, energy and capacity based incentives for carbon-free resources, renewable portfolio standards (RPS), and many others \cite{hogan2015cleaner}. A wide literature has shown that carbon pricing is the optimal tool for efficiently decreasing $\mathrm{CO_2}$ emissions from the electricity sector as well as from the whole economy \cite{aldy2010designing,fischer2008environmental,palmer2011federal}. Ideally, carbon pricing would be the main decarbonization instrument, with $\mathrm{CO_2}$ emissions priced according to the marginal damages and low-emission resources (e.g. renewables, storage, demand-side management) given necessary subsidies. For example, these subsidies can be based on the R\&D spillover rate or compensate for generation spillovers from learning-by-doing \cite{fischer2008environmental}. However, putting a price on a unit of $\mathrm{CO_2}$ emission has been politically unpopular\footnote{In the U.S., carbon pricing initiatives have failed both on the federal level (e.g. Waxman-Markey bill, Coons/Feinstein bill and many others) and state level (I-732 initiative and SB 5971 bill in Washington, SB1530 bill in Oregon, among others). In France, carbon pricing legislation was one of the factors behind the popular protest movement of the “Yellow Vests.” See also \cite{Soren2019} for an attempt to explain why carbon taxes are hard to pass.} such that in practice, only few jurisdictions maintain meaningful carbon prices \cite{Dsire}, i.e. put a price on $\mathrm{CO_2}$ emissions that comes close to the value of marginal damages from these emissions. Instead, policy-makers have widely chosen to rely on second-best policies of subsidizing carbon-free generation and investment (see \cite{hogan2015cleaner} for some  of the unintended consequences of subsidies).

Ill-designed policies can distort electricity markets and  benefit certain stakeholders over others,  even leading  to artificially prolonged economic life of certain generation types, and   result in inefficiencies of the entire market and power system. Entangled interests of stakeholders and their strategic behavior in the electricity market, along with the physical interconnection through the system, complicate the decision-making process even further. Hence,  there is a need in policy design tools, which are capable of taking into account different techno-economic perspectives of stakeholders   (e.g. regulators, power utilities, generation companies, consumers, etc.)   and physical limits on power system operations.

Among a diverse set of policy instruments for grid decarbonization \cite{hogan2015cleaner}, we narrow our attention to incentives for carbon-free resources, such as RES, given the relative prevalence and success of these incentives in the US and globally. The RES incentives fall into two categories: energy- and capacity-based. The former policies remunerate RES producers proportionally to their actual renewable production (e.g. FIT, REC), while the latter policies remunerate RES producers based on their installed capacity (e.g. PTC). Many countries use either   energy or capacity-based policies, or both of them, depending on their socio-economic circumstances \cite{murdock2019renewables}. The majority of U.S. states support RPS goals using RECs or PTCs, while Great Britain, Japan, and China have recently turned to FITs \cite{murdock2019renewables}. Regardless of the type, these policies have played  a leading role in rolling out 1,449 GW of renewables all over the world (260 GW in the U.S.) \cite{ajadi2019global}.

Methodologically, the literature on RES incentives has  been dominated by empirical data analyses, rather than model-based analyses. Huntington et al. \cite{huntington2017revisiting} compared the performance of capacity- and energy-based incentives in Spain and concluded that capacity-based incentives are more compatible with electricity markets in the considered case because they do not distort  market outcomes and generation payments, thus  not affecting the  efficiency of market competition in the considered setting. Furthermore, Newbery et al.\cite{newbery2018market} raised concerns that FITs, which are a non-competitive revenue stream that drives  RES profitability,  cause locational  market distortions affecting conventional producers. On the other hand, \"Ozdemir et al.\cite{ozdemir2019capacity} demonstrated that energy-based incentives are more cost-effective for incentive providers to achieve production goals, whereas  capacity-based incentives are more desirable in the longer term and allow for achieving  RES-rich RPS goals. Furthermore, Nicolini and Tavoni \cite{nicolini2017renewable} analyzed the effectiveness of RES incentives in five European countries and, in line with  \cite{ozdemir2019capacity}, concluded that FITs outperformed tradable green certificates (similar to RECs in the U.S.) in terms of achieving a greater RES production over the period from 2000 to 2010. Although both capacity- and energy-based RES incentives have pros and cons that have been demonstrated via empirical data analyses in \cite{huntington2017revisiting, newbery2018market, ozdemir2019capacity, nicolini2017renewable},  regulatory and legislative bodies still lack decision-support tools for (i) exogenous policymaking (i.e., optimization of incentives) and (ii) on-par comparison of different incentive types.

A critical point in designing renewable energy incentives is the need to accommodate and balance the perspectives of multiple stakeholders, as for example state regulators, power utilities, generation companies, consumers, etc. Ideally, regulators aim to select such electricity tariffs and policy incentives that maximize the social welfare in their jurisdiction, e.g. to  achieve a given RPS target in the least-cost manner, and power utilities seek to profit from supplying electricity to consumers, while ensuring  reliable electricity supply. Naturally, these objectives are conflicting because  power utilities  benefit from a higher electricity tariff, while  a higher electricity tariff leads to a greater consumer payment, which in turn reduces the social welfare. Furthermore, the wholesale electricity market typically solves a social-welfare maximization or cost minimization problem, which again trade offs the  interests of power utilities and generation companies. Thus, setting renewable energy incentives without accounting for the conflicting interests of stakeholders may yield both  market distortions and inefficient renewable expansion plans. This motivates the need for developing a decision-support tool to optimally set the electricity tariff and incentives, which balances the  interests of each stakeholder.

This paper takes the perspective of the state regulator and develops a modeling framework to jointly optimize RES incentives and electricity tariffs, {under the given RPS goal and incentive budget}, while incorporating interactions between the stakeholders involved. To appropriately model this multi-stakeholder environment, the proposed framework is formulated as a tri-level (TL) optimization problem, where each level represents a different entity: a state regulator, a power utility, and a wholesale electricity market. Additionally, the level representing the  power utility internalizes flexible demands and  time-of-use (TOU) tariffs to assess the impact and value of the demand-side flexibility on achieving the RPS goals. Using this TL optimization, this paper studies a complex trade-off between the RES incentives and electricity tariffs under RPS goals and different retirement scenarios.

In order to realistically model power system operations and, in particular, the need in balancing capacity to support the roll-out of high penetration RES levels, we model the RES stochasticity using chance constraints, as in \cite{bienstock2014chance}, which is shown to have superior computational performance over traditional scenario-based stochastic programming. However, even if chance constraints are used to alleviate some computational complexity,  commercial solvers are still not able to deal with  multi-level structure problems efficiently and, therefore, algorithmic adaptations are inevitable. Thus, to solve the proposed TL problem, we use the column-and-cut generation (C\&CG) method,  which decomposes the original problem into master and sub-problems and solves them iteratively until they reach convergence. Unlike previous C\&CG implementations in \cite{ruiz2015robust,duarte2020multi,dvorkin2017co,moreira2014adjustable}, the proposed TL problem  contains bilinear terms with the multiplication of multiple-level variables. Thus, we enhance the original C\&CG algorithm further to contain auxiliary columns for generating cuts accounting for bilinear terms.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.82\columnwidth]{Figures/ModelStructure_updated.pdf}
    \caption{\small Structure of the proposed TL model. Note that $\ast$ denotes a parameterized value of decision variables from lower decision-making levels.}
    \label{fig:model}
\end{figure}

\section{Model}\label{Sec:model}

Resorting to the TL structure, Fig.~\ref{fig:model} illustrates the proposed  model with three different stakeholders: a state regulator, a power utility and a wholesale electricity market. Each entity is represented by an optimization problem and the TL framework makes it possible to implement a leader-follower structure among entities, i.e. the state regulator is strategically more dominant than the power utility, and the power utility can act strategically in the wholesale electricity market. The goal of the upper-level (UL) problem is to represent a decision-making process of the state regulator, which sets the electricity tariffs and RES incentives and monitors performance and policy compliance of power utilities.  The middle-level (ML) problem parameterized with the UL decisions represents a decision-making process of a typical US electric power distribution utility, which decides  which generations assets, both RES and conventional, to build in order to supply electricity to its consumers and how to operate newly built and existing generation assets given representative load and renewable forecasts. The ML operational decisions include the dispatch of utility-owned assets and its participation strategy in the wholesale electricity market. Each state is assumed to have one power utility connected to the transmission network via the  power flow interface with the  limit of $P_n^{\downarrow, \max}$. Finally, the lower-level (LL) problem represents the wholesale electricity market, in which the power utility modeled in the ML problem acts strategically and the generation companies and power utilities from other states are  price takers.

\subsection{UL problem (State regulator)}

In the U.S., states can set renewable policies and provide incentives ($\tau^{\mathrm e}, \tau^{\mathrm c}$) to support their implementation. At the same time, they regulate electricity tariffs ($\pi_{nt}$) so that the consumer payments are not unjustly manipulated by power utilities, while power utilities are guaranteed to recover the incurred costs to supply electricity reliably.
Following this practice, we formulate the problem of the state regulator as:
\begin{subequations}
\begin{align}
    \begin{split}
        & \min_{\Xi^{\mathrm{R}}}  O^{\mathrm{R}}_s \coloneqq \hspace{-0.5mm} \sum_{t\in{\cal T}}\sum_{e\in{\cal E}} \omega_e \Big[ \sum_{n\in{\cal N}_s} D_{nte} \pi_{nt}
        + \hspace{-2.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}\cup{\cal I}^{\mathrm R}} \hspace{-3mm}g_{ite} \tau^{\mathrm{e}}  \Big]
        \\&\hspace{20mm}
        + \hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}}  g^{\mathrm{max}}_i  \tau^{\mathrm{c}}
    \end{split}\label{Eq:RegObj}\\
     & \sum_{t \in {\cal T}} \sum_{i\in{\cal I}_s^{\mathrm R}\cup\hat{\cal I}_s^{\mathrm R}} \!\!\!\!{g}_{ite} \tau^{\mathrm{e}} +   \sum_{i\in\hat{\cal I}_s^{\mathrm R}} {g}^{\mathrm{max}}_{i} \tau^{\mathrm{c}}  \le  B^{\mathrm{P}}_s,~\forall e\in{\cal E},\label{Eq:Tariff_budget}\\
     & \sum_{t \in {\cal T}} \sum_{i \in {\cal I}_s^{\mathrm R}\cup \hat{\cal I}_s^{\mathrm R}} \!\!{g}_{ite} \ge \kappa_s \sum_{t \in {\cal T}} \sum_{n\in{\cal N}_s} {D}_{nte},~\forall e\in{\cal E}, \label{Eq:StrReg_kappa}\\
     & O^{\mathrm{U}}({g}^{\max}_i, {g}_{ite}) \ge 0,\label{Eq:Radeq}\\
     & \pi_{nt} =
     \begin{cases}
        &\pi^{\mathrm{on}}_n,\quad\forall t \in {\cal T}^{\mathrm{on}},\\
        &\pi^{\mathrm{off}}_n,\quad\forall t \in {\cal T}^{\mathrm{off}},
    \end{cases}\label{Eq:TOU_model1}
\end{align}\label{ULproblem}\end{subequations}
where $\Xi^{\mathrm{R}}=\{\pi_{nt}, \pi^{\mathrm{off}}_n, \pi^{\mathrm{on}}_n, \tau^{\mathrm{e}}, \tau^{\mathrm{c}} \ge 0\}$. Equation~\eqref{Eq:RegObj} minimizes the expected total consumer payments and the total incentives provided to support renewable expansion over representative days $e$ with probability $\omega_e$. This formulation is based on two assumptions on the state regulator. First, the state regulator aims to meet an explicitly given RPS goal, while keeping the  total cost  incurred by power sector as low as possible. Second, the  state regulator can raise money for subsidies without creating welfare losses, i.e. that the marginal cost of public funds equals zero. Also,  this formulation considers  a particular case of a multi-objective problem with the equal weight between the consumer payment and the total incentives, whereas this weighting preference can be tuned according to the preference of the state regulator  and according to the actual costs of raising public funds. The capacity-based incentive ($\tau^{\mathrm c}$) has a net present value  prorated on a daily basis. The renewable policy budget limit is imposed in Eq.~\eqref{Eq:Tariff_budget}, and is procured externally from the optimization and is decoupled from the electricity tariff ($\pi_{nt}$). The RPS target constraint, representing a given renewable policy, is in Eq.~\eqref{Eq:StrReg_kappa} and the revenue adequacy (non-negative profit) constraint for the power utility is imposed in Eq.~\eqref{Eq:Radeq} to ensure the power utility is not exposed to the deficit by investing in RES, where $O^{\mathrm{U}}$ is the objective function of the utility defined in Eq.~\eqref{Utility_Obj}. The non-negativity constraint on profits of the utility does not preclude retirement of the assets -- as we do not model decommission costs or annual, fixed maintenance costs, retirement of a power plant is equivalent to having zero capacity utilization rate of that power plant. The electricity tariff is modeled as a time-of-use rate in Eq.~\eqref{Eq:TOU_model1} while ${\cal T}^{\mathrm{on/off}}$ represent on- and off-peak time sets.

\subsection{ML problem (Power utility)}
The ML problem models the power utility.  First, we consider  flexible demand and  generation models.

\subsubsection{Demand} Demand elasticity  can be implicitly modeled in the ML optimization by deriving optimality conditions of the surplus maximization problem of aggregated consumers. An aggregated utility function of electricity consumers is  non-decreasing and concave and modeled as \cite{samadi2010optimal}:
\begin{subequations} \begin{align}
    & U(d_{nte}) =
        \begin{cases}
            M d_{nte} - \frac{N}{2} d_{nte}^2 \hspace{3mm} \text{if} \quad 0\le d_{nte} \le \frac{M}{N},\\
            \frac{M}{N} \hspace{22.7mm} \text{if} \quad  d_{nte} \ge \frac{M}{N},
        \end{cases}
\end{align}\end{subequations}
where $M$ and $N$ are pre-determined parameters.
Given electricity tariff $\pi_{nt}$, the maximized consumer surplus is:
\begin{subequations} \label{eq:flexible_demand}\begin{align}
    & \max\hspace{1mm} \text{CS}_{nte} \coloneqq(M d_{nte} - \frac{1}{2} N d_{nte}^2) - \pi_{nt} d_{nte}
\end{align}
Then, we derive a closed-form demand  as a function of the electricity tariff ($d_{nte}(\pi_{nt})$) via the first-order optimality condition and by assuming non-zero marginal utility of electricity (setting upper bound $\overline{d}_{nte}=M/N$):
\begin{align}
    &\frac{\partial {\text{CS}_{nte}}}{\partial d_{nte}} = (M - \pi_{nt}) -N d_{nte} = 0,
    \\&d_{nte}(\pi_{nt}) = \frac{M - \pi_{nt} }{N}.  \label{eq:demand_simple}
\end{align}\end{subequations}
Given \eqref{eq:demand_simple},  the demand function depends on parameters $M$ and $N$, which are external to the optimization problem, and the value of tariff chosen by the UL optimization. If the power utility implements a demand-response program with incentives for consumers to change their power consumption, the resulting demand in \eqref{eq:demand_simple}  can be further adjusted. For instance, if $\Delta \pi_{nt} \geq 0$ denotes the demand-response incentive, the resulting demand in \eqref{eq:demand_simple}  will be set to $d_{nte}(\pi_{nt} + \Delta \pi_{nt}) = \frac{M - \pi_{nt} - \Delta \pi_{nt} }{N}.$

\subsubsection{Generation}
The uncertain RES generation output is modeled through forecast error variables $\bm{\epsilon}_{ite}$ with Gaussian distribution \footnote{The Gaussian distribution is used in this work due to its mathematical tractability, while other probability distributions  can be adopted to model  forecast errors with more accuracy (e.g. Laplace, logistic, student's t-distributions, see \cite{roald2015security}).}
${\cal G}({\mathbb{E}(\bm{\epsilon}_{ite})}, {\mathrm{Stdev}(\bm{\epsilon}_{ite})}^2)$ as:
\begin{subequations}\label{Eq:Random_model}\begin{align}
\bm{\epsilon}_{ite} \!\!\sim\!\! \begin{dcases}
    &\hspace{-3mm} {\cal G}({G^{\mathrm{max}}_i}\upsilon_{ite}, ({G^{\mathrm{max}}_i} \sigma_{ite})^2) ,\hspace{0.0mm}\forall i \in{\cal I}^{\mathrm R}_s, t \in {\cal T}\!, e\in{\cal E},\hspace{-3mm}
    \\&\hspace{-3mm} {\cal G}({g^{\mathrm{max}}_i}\upsilon_{ite}, ({g^{\mathrm{max}}_i} \sigma_{ite})^2),\hspace{2mm}\forall i \in\hat{\cal I}^{\mathrm R}_s, t \in {\cal T}\!, e\in{\cal E},\hspace{-3mm}
\end{dcases}
\end{align}
where $\upsilon_{ite}$ and $\sigma_{ite}$ are the normalized mean and standard deviation of the forecast error, and ${G^{\mathrm{max}}_i}$ and ${g^{\mathrm{max}}_i}$ are the generation capacity of existing (${\cal I}^{\mathrm R}_s$) and candidate ($\hat{\cal I}^{\mathrm R}_s$) RES generators.
Next, the forecast power output of RES generators, $\overline{g}_{ite}$, can be modeled with a forecast factor $\rho_{ite}\in [0,1]$ as:
\begin{align}
\overline{g}_{ite} \!=\! \begin{dcases}
    &\hspace{-3mm} \rho_{ite} G_i^{\mathrm{max}} ,\hspace{5mm}\forall i \in{\cal I}^{\mathrm R}_s,t \in {\cal T}, e\in{\cal E},
    \\&\hspace{-3mm} \rho_{ite} g_i^{\mathrm{max}},\hspace{6mm}\forall i \in\hat{\cal I}^{\mathrm R}_s, t \in {\cal T}, e\in{\cal E},
\end{dcases}
\end{align}
Then, the output of RES generators follows as:
\begin{align}
    & \bm{g}_{ite} = \overline{g}_{ite} + \bm{\epsilon}_{ite}, \quad\forall i \in {\cal I}^{\mathrm{R}}_s\cup\hat{\cal I}^{\mathrm{R}}_s, t \in {\cal T}, e\in{\cal E}, \label{Eq:uncertainty_model_n}\vspace{-3mm}
\end{align}
Since forecast errors result in a real-time mismatch between the power produced and consumed, controllable generation resources must offset these deviations. In practice, the real-time affine control is used, \cite{jaleeli1992understanding}, and modeled as:
\begin{align}
    & \bm{g}_{ite} = \overline{g}_{ite} \!-\! \alpha_{it} \!\!\!\!\!\!\!\!\!\!\sum_{j \in {\cal I}^{\mathrm R}_{s(i)}\cup\hat{\cal I}^{\mathrm R}_{s(i)}} \!\!\!\!\!\!\!\!\!\!\bm{\epsilon}_{jte},~\forall i \in {\cal I}^{\mathrm C}_s\!\cup\!\hat{\cal I}^{\mathrm C}_s, t \in {\cal T}, e\in{\cal E}, \label{Eq:affine_control}
\end{align}%
\end{subequations}%
where $\sum_{i\in{\cal I}^{\mathrm C}_s\cup\hat{\cal I}^{\mathrm C}_s} \alpha_{it} \!=\! 1$ ensures the sufficiency of procured balancing resources. Synchronized generators continuously offset the power mismatch based on  parameter $\alpha_{it}$, which regulates the participation of each generator in balancing. The value of $\alpha_{it}$ can be     set ahead of time, e.g. $\alpha_{it}=1/{\text{card}({\cal{I}}^{\mathrm{C}}_s\cup \hat{\cal{I}}^{\mathrm{C}}_s})$, or optimized as in \cite{bienstock2014chance}.

\subsubsection{Power utility model}
We assume the power utility company optimizes its generation expansion plan and optimizes its operations, including the participation strategy in the wholesale electricity market as follows:
\begin{subequations}\begin{align}
    \begin{split}
        &\max_{\Xi^{\mathrm{U}}} O^{\mathrm{U}}\!\coloneqq\!
        \hspace{-0.0mm}\sum_{e\in{\cal E}} \omega_e \Big[\sum_{t\in{\cal T}}\sum_{n\in{\cal N}_s}\pi_{nt}d_{nte}
        + \sum_{t\in{\cal T}}\hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}\cup{\cal I}^{\mathrm R}}  \hspace{-3mm}\tau^{\mathrm{e}} g_{ite}   \hspace{-5mm}
        \\&\hspace{0mm}
        + \hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}} \tau^{\mathrm{c}} g^{\mathrm{max}}_i
        - \hspace{-1mm}\sum_{t\in{\cal T}} \Big( \hspace{-0.5mm}\sum_{n\in{\cal N}_s} \hspace{-1.5mm} \lambda_{nte} d_{nte}\hspace{-1mm}
        -  \hspace{-2mm}\sum_{i \in {\cal I}\cup\hat{\cal I}} \hspace{-1.5mm}\lambda_{n(i),t,e} {g}_{ite} \!\Big) \hspace{-5mm}
        \\&\hspace{0mm}
        - \sum_{i\in\hat{\cal I}} C^{\mathrm{inv}}_i g^{\mathrm{max}}_{i}
        - \sum_{i\in{\cal I}}\sum_{t\in{\cal T}} C^{\mathrm{g}}_i g_{ite}\Big]\hspace{-5mm}
        \label{Utility_Obj}
    \end{split}\\
    & \mathbb{P}[\bm{g}_{ite} \le {G}_{i}^{\mathrm{max}}]\ge1-\eta,\quad\forall i \in {\cal I}^{\mathrm{C}}_s, t \in {\cal T}, e\in{\cal E},\label{Eq:gen_UB}\\
    & \mathbb{P}[\bm{g}_{ite} \ge {G}_{i}^{\mathrm{min}}]\ge1-\eta,\quad\forall i \in {\cal I}^{\mathrm{C}}_s, t \in {\cal T}, e\in{\cal E}, \label{Eq:gen_LB}\\
    & \mathbb{P} [\bm{g}_{ite} \le {g}_{i}^{\mathrm{max}}] \ge 1 - \eta,\quad\forall i \in {\cal \hat{I}}^{\mathrm{C}}_s, t \in {\cal T}, e\in{\cal E},  \label{Eq:gen_UB_new}\\
    & \mathbb{P} [ \bm{g}_{ite} \ge \Gamma_i {g}_{i}^{\mathrm{max}}] \ge 1 - \eta,\quad\forall i \in {\cal \hat{I}}^{\mathrm{C}}_s, t \in {\cal T}, e\in{\cal E},  \label{Eq:gen_LB_new}\\
    & {H}^{\mathrm{min}}_{i} \!\!\le\! \overline{g}_{ite}\!-\!\overline{g}_{i,t-1,e} \! \le \!{H}^{\mathrm{max}}_{i}\!\!,~\!\forall i \!\in \! {\cal I}^{\mathrm{C}}_s\!\cup\!\hat{\cal I}^{\mathrm{C}}_s, t \in \!{\cal T}\!, e\in\!{\cal E},\hspace{-3mm}\label{Eq:ramp}\\
    & \sum_{i \in {\cal I}_s \cup {\hat{\cal{I}}_s}} \!\!\!\overline{g}_{ite} + \!{p}_{nte}^{\downarrow}  = \!d_{nte} ,\quad\forall t \in {\cal T},{n\in{\cal N}_s}, e\in{\cal E}, \label{TSO_P_balance2}\\
    & -\!{P}_n^{\downarrow,\mathrm{max}} \!\!\le {p}^{\downarrow}_{nte} \le {P}_n^{\downarrow,\mathrm{max}},~ \forall t \in {\cal T}\!, n\in{\cal N}_s, e\in{\cal E},\label{Eq:pdown_ULB}
\end{align}\label{MLproblem}\end{subequations}
where ${\Xi^{\mathrm{U}}=\{g^{\mathrm{max}}_i, \overline{g}_{ite} \ge 0\}}$. Objective function \eqref{Utility_Obj} maximizes the  profit of the power utility over the operating horizon $t\in{\cal T}$ for each representative day $e$, which includes the payment collected from consumers, revenue from the wholesale market participation and renewable energy incentives provided by the state regulator. Investment cost $C^{\mathrm{inv}}_i$ and capacity-based incentive $\tau^{\mathrm c}$ are prorated on a daily basis using the net present value approach as in \cite{6936940}. Eqs.~\eqref{Eq:gen_UB}-\eqref{Eq:gen_LB_new} limit power outputs of existing and candidate resources, while the RES uncertainty is captured via chance constraints with random variable $\bm{g}_{ite}$. Eq.~\eqref{Eq:ramp} limits the up- and downward ramping rates for controllable generators. Finally, Eq.~\eqref{TSO_P_balance2} enforces the supply demand balance, while Eq.~\eqref{Eq:pdown_ULB} enforces the power flow limit ($p^{\downarrow}_{nte}$) on the power flow between the power utility and transmission network (and  the  wholesale electricity market). Note that network constraints among nodes $n\in{\cal N}_s$ within each state $s$ are not imposed in Eq.~\eqref{MLproblem}, and the power flows are constrained in the wholesale market problem, see Eq.~\eqref{LLproblem}.

\subsection{LL problem (Wholesale market)}
The wholesale electricity market for each representative day $e\in{\cal E}$, which includes the strategic participant  modeled in the ML problem and other non-strategic participants, solves:
\begin{subequations}\begin{align}
    \begin{split}
        &\max_{\Xi^{\mathrm{W}}} O^{\mathrm{W}}_{e}\!\coloneqq\!
       \sum_{t\in{\cal T}} \sum_{i \in {\cal I}\cup\hat{\cal I}} - C_i^{\mathrm{g}} {g}_{ite}\label{WM_Obj}
    \end{split}\\
    &(\xi_{lte}):\!\!\!\quad f_{lte} = \frac{1}{X_l} (\theta_{o(l),te}-\theta_{r(l),te}),~\forall {l \in {\cal L}},~t\in{\cal T},\label{TSO_DSPF}\\
    \begin{split}
        &(\lambda_{nte}):\hspace{-1mm}\sum_{i \in {\cal I}_n\cup\hat{\cal I}_n } \hspace{-1mm}{g}_{ite}  + \hspace{-1mm}\sum_{l \vert r(l) = n } \hspace{-1mm} f_{lte} \!-\hspace{-1mm}\sum_{l \vert o(l) = n } \hspace{-1mm} f_{lte} \!
        \\&\hspace{26mm}=
        \begin{dcases}
            &\hspace{-3mm}d_{nte},\quad\forall n\in{\cal N}_s,~t\in{\cal T},\\
            &\hspace{-3mm}D_{nte},\quad\forall n\in{\cal N}\backslash{\cal N}_s,~t\in{\cal T},
        \end{dcases}
    \end{split}\label{TSO_P_balance}\\
    &(\underline{\gamma}_{ite},\overline{\gamma}_{ite}):\!\quad 0 \le {g}_{ite} \le \overline{g}_{ite},~\forall i \in {\cal{I}}\cup\hat{\cal{I}},~t\in{\cal T}\!,\label{TSO_G_plim}\\
    &(\underline{\delta}_{lte},\overline{\delta}_{lte}):\!\quad -{F}^{\mathrm{max}}_l \le f_{lte} \le {F}^{\mathrm{max}}_l,~\forall {l \in {\cal L}},~t\in{\cal T}\!,\!\!\label{TSO_LineLim}
\end{align}\label{LLproblem}\end{subequations}
where~$\Xi^{\mathrm{W}}=\{g_{ite}, \theta_{nte}, f_{lte} \ge 0\}$. Objective function \eqref{WM_Obj} maximizes the social welfare under inflexible wholesale demand\footnote{Note that the flexible demand modeled in Eq.~\eqref{eq:flexible_demand} is internalized in the ML problem, i.e. the power utility leverages all flexibility, and the wholesale market thus deals with fixed demand curves, denoted as $d_{nte}$ for nodes in strategically acting state and $D_{nte}$ for others.}, i.e. it minimizes the operating cost. Equation~\eqref{TSO_DSPF} models DC power flows and Eq.~\eqref{TSO_P_balance} enforces the nodal power balance in the transmission system, where $o(l)$ and $r(l)$ denote originating and receiving nodes of line $l$, respectively. The dual variable of Eq.~\eqref{TSO_P_balance} ($\lambda_{nte}$) is then the locational marginal price (LMP). Equations~\eqref{TSO_G_plim}-\eqref{TSO_LineLim} limit generation outputs and transmission line flows.

\section{Solution Method}
To solve the proposed TL model effectively, we convert chance constraints into a deterministic form, derive an equivalent bi-level problem, and develop a C\&CG algorithm variant that deals with bilinear terms across levels.

\subsection{Chance constraint reformulation}\label{subsec:CCreform}
As in \cite{kim2020computing}, chance constraints in Eqs.~\eqref{Eq:gen_UB}-\eqref{Eq:gen_LB_new} can be exactly reformulated as the following deterministic constraints:
\begin{subequations}\label{Eq:CC_reform}
\begin{align}
    \begin{split}
        & \overline{g}_{ite}
        -\alpha_{it} \Big( \hspace{-0.0mm} \sum_{j\in{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} G^{\mathrm{max}}_j \upsilon_{jte}
        \!+ \hspace{-1.0mm} \sum_{j\in\hat{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} g^{\mathrm{max}}_j \upsilon_{jte}\Big)
        \\& \hspace{5.5mm}+ \Phi^{-1}(1\!-\!\eta){\Stdev[\bm{g}_{ite}]}
        \le {G}_{i}^{\mathrm{max}},~\forall i \in  {\cal I}^{\mathrm{C}}_s, t \in {\cal T},\!\!
    \end{split}\label{Eq:gen_UB_Conic}\\
    \begin{split}
        & \overline{g}_{ite}
        -\alpha_{it} \Big( \hspace{-0.0mm} \sum_{j\in{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} G^{\mathrm{max}}_j \upsilon_{jte}
        \!+ \hspace{-1.0mm} \sum_{j\in\hat{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} g^{\mathrm{max}}_j \upsilon_{jte}\Big)
        \\& \hspace{5.5mm}+ \Phi^{-1}(1\!-\!\eta){\Stdev[\bm{g}_{ite}]}
        \le {g}_{i}^{\mathrm{max}},~\forall i \in  \hat{\cal I}^{\mathrm{C}}_s, t \in {\cal T},
    \end{split}\label{Eq:gen_UB_New_Conic}\\
    \begin{split}
        & \overline{g}_{ite}
        -\alpha_{it} \Big( \hspace{-0.0mm} \sum_{j\in{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} G^{\mathrm{max}}_j \upsilon_{jte}
        \!+ \hspace{-1.0mm} \sum_{j\in\hat{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} g^{\mathrm{max}}_j \upsilon_{jte}\Big)
        \\& \hspace{5.5mm}- \Phi^{-1}(1\!-\!\eta){\Stdev[\bm{g}_{ite}]}
        \ge {G}_{i}^{\mathrm{min}},~\forall i \in  {\cal I}^{\mathrm{C}}_s, t \in {\cal T}, \!
    \end{split}\label{Eq:gen_LB_Conic}\\
\end{align}
\begin{align}
    \begin{split}
        & \overline{g}_{ite}
        -\alpha_{it} \Big( \hspace{-0.0mm} \sum_{j\in{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} G^{\mathrm{max}}_j \upsilon_{jte}
        \!+ \hspace{-1.0mm} \sum_{j\in\hat{\cal I}^{\mathrm{R}}_{s}} \hspace{-1.0mm} g^{\mathrm{max}}_j \upsilon_{jte}\Big)
        \\& \hspace{5.5mm}- \Phi^{-1}(1\!-\!\eta){\Stdev[\bm{g}_{ite}]}
        \ge \Gamma_i {g}_{i}^{\mathrm{max}},~\forall i \in  \hat{\cal I}^{\mathrm{C}}_s, t \in {\cal T}, \hspace{-10mm}
    \end{split}\label{Eq:gen_LB_New_Conic}\\
    & \Stdev[\bm{g}_{ite}] \!=\! \alpha_{it}\!\sqrt{\sum_{j\in{\cal{I}}_{s(i)}^{\mathrm{R}}}({{G}^{\mathrm{max}}_{j}} \sigma_{jte})^2 \!+\!\! \sum_{j\in\hat{\cal{I}}_{s(i)}^{\mathrm{R}}}({{g}^{\mathrm{max}}_{j}} \sigma_{jte})^2 }, \label{Eq:stdev}
\end{align}\end{subequations}
where $\Phi^{-1}(\cdot)$ is the inverse cumulative distribution function of the Gaussian distribution and $\Stdev[\cdot]$ is a standard deviation operator. Eq.~\eqref{Eq:stdev} assumes independence of random variables $\{\bm{g}_{ite} ~| ~\forall i\!\in \!{\cal I}_s^{\mathrm{R}}\!\cup\!\hat{\cal I}_s^{\mathrm{R}}\}$ across time $t$ and representative days $e$.

\subsection{Equivalent bi-level problem}\label{subsec:bilevelreform}
Since the LL problem in Eq.~\eqref{LLproblem} is a linear problem, it is guaranteed to have a strict complementary solution by the \textit{Goldman-Tucker theorem} \cite{goldman1956theory}. Hence, we can apply the strong duality theorem to the LL problem and  convert it into a set of constraints, resulting in an equivalent ML and LL problem. This equivalent problem makes it possible to convert the proposed TL problem into an equivalent bi-level problem.

First, we formulate the dual problem of Eq.~\eqref{LLproblem} for each representative day $e\in{\cal E}$ as:

\begin{subequations}\begin{align}
    \begin{split}
        & \!\underset{\Xi^{\mathrm{DW}}}{\min} O^{\mathrm{DW}}_{e}\coloneqq \sum_{t\in {\cal T}} \Big( \sum_{i \in {\cal I}\cup\hat{\cal I}} \! \overline{g}_{ite} \overline{\gamma}_{ite}
        - \!\!\sum_{n \in {\cal N}} \!D_{nte}^{\mathrm{p}} \lambda_{nte} \!\!\\
        &\hspace{30mm}   + \sum_{l \in {\cal L}} \! \big({F}^{\mathrm{max}}_l \overline{\delta}_{lte}\!+\!{F}^{\mathrm{max}}_l \underline{\delta}_{lte} \big) \Big)
    \end{split}\label{TSO_Obj_dual}\\
    & (g_{ite}^{\mathrm p}):\hspace{2mm}-\underline{\gamma}_{ite}\!+\!\overline{\gamma}_{ite}\!-\!\lambda_{n(i),te} \!=\! - C^{\mathrm g}_i ,~\forall i \in {\cal{I}},t \in {\cal T},\label{DLL_gp}\\
    \begin{split}
        & (f_{lte}^{\mathrm p}):\hspace{2mm}\!-\underline{\delta}_{lte}+\overline{\delta}_{lte}+\xi_{lte} \!-\!\lambda_{r(l),t,e}+\lambda_{o(l),t,e} = 0 ,
        \\&\hspace{55mm}\forall l \in {\cal{L}},t \in {\cal T},
    \end{split}\label{DLL_fl}\hspace{-1mm}\\
    & (\theta_{nte}):\hspace{2mm}-\!\!\!\!\sum_{l \vert o(l) = n}\! \frac{\xi_{lte}}{X_l}\!+\!\!\!\!\!\sum_{l \vert r(l) = n}\! \frac{\xi_{lte}}{X_l} = 0, ~\forall n \in {\cal{N}},t \in {\cal T},\label{DLL_theta}\!
    \end{align}\end{subequations}
where $\Xi^{\mathrm{DW}} = \big\{ \underline{\gamma}_{ite}, \overline{\gamma}_{ite}, \underline{\delta}_{lte}, \overline{\delta}_{lte} \geq 0;  \xi_{lte}, \lambda_{nte} \text{: free} \big\}$.

Then, the proposed TL problem can be equivalently recast as a compact bi-level problem:
\begin{subequations}\label{Eq:bilevel}\begin{align}
    & \text{Eq.}~\eqref{ULproblem}\hspace{31.5mm}:\text{UL problem}\\
    \begin{split}
        & \Xi^{\mathrm{U}} \cup \Xi^{\mathrm{W}} \in \arg\big[
        \\& \hspace{7mm}\text{Eqs.}~\eqref{Utility_Obj},\eqref{Eq:ramp}\text{\textendash}\eqref{Eq:pdown_ULB}, \eqref{Eq:CC_reform}\hspace{1.5mm} :\text{ML problem}\\
        & \hspace{7mm}\text{Eqs.}~\eqref{TSO_DSPF}\text{\textendash}\eqref{TSO_LineLim}\hspace{14mm}:\text{LL primal constraints}\\
        & \hspace{7mm}\text{Eqs.}~\eqref{DLL_gp}\text{\textendash}\eqref{DLL_theta}\hspace{14mm}:\text{LL dual constraints}\\
        & \hspace{7mm}O^{\mathrm{W}}_{e}=O^{\mathrm{DW}}_{e},\quad\forall e\in{\cal E}\hspace{1.3mm}:\text{LL strong duality}
        \hspace{3mm}\big].
    \end{split}\label{Eq:combinedMLLL}
\end{align}\end{subequations}
The complete formulation of \eqref{Eq:bilevel} with all constraints detailed is described in Appendix.
It is noteworthy that duality theory is not applicable to the combined ML and LL problem in \eqref{Eq:combinedMLLL} due to the existence of bilinear terms $\tau^{\mathrm{e}} g_{ite},\tau^{\mathrm{c}} g_{i}^{\mathrm{max}}, \lambda_{n(i),t,e} g_{ite}$ in objective functions, and therefore, further conversion of the bi-level problem in \eqref{Eq:bilevel} into an exact single-level form is unattainable. This motivates the use of the C\&CG algorithm for solving \eqref{Eq:bilevel} described in the next subsection.

\subsection{Column-and-cut generation algorithm}
\begin{figure}[b]
    \centering
        \includegraphics[width=0.9\columnwidth]{Figures/CCG_SD.pdf}
        \caption{\small Implementation of the  C\&CG algorithm. \vspace{-0mm}}
        \label{fig:CCGstructure}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.98\columnwidth]{Figures/CCG_alg_updated.pdf}
    \vspace{-0mm}
\end{figure}

The C\&CG algorithm exploits the primary-secondary (also known as master-slave) structure of Eq.~\eqref{Eq:bilevel} and generates primal cuts and columns over iterations using primal information exchanges between the master- and sub-problems. In our case, as shown in Fig.~\ref{fig:CCGstructure}, the master problem provides policy decisions of the state regulator to the subproblem and receives reported performance (expansion and dispatch decisions) in return.

The C\&CG algorithm as applied to Eq.~\eqref{Eq:bilevel} is detailed in Algorithm 1. Step 1 initializes the lower and upper bounds of the objective function ($LB,UB$) and the iteration counter ($k$). Step 2-1 obtains initial solutions of the master problem without any cut and update $LB$ and given the UL solution, while Step 2-2 solves a subproblem, saves solutions, and update $UB$. Step 3 first updates the iteration counter (Step 3-1). Then  it obtains the master problem solution  with $k\!-\!1$ sets of cuts generated from the subproblem solutions from all the previous iterations ($1,\dots,k\!-\!1$) and updates $LB$ (Step 3-2). Finally, it solves the subproblem with the solution of the master problem at the current iteration $k$ and updates $UB$ (Step 3-3). This process repeats until the convergence is achieved at a desired tolerance level, $(UB-LB)/UB \le \epsilon$, i.e., the approximation error of the master problem on the subproblem objective $(O^{\mathrm U})$ is within the desired tolerance level.

\subsubsection{Subproblem}
Equation~\eqref{SP} solves the combined ML and LL equivalent problem, while the UL decisions of the current iteration $\Xi^{\mathrm{R}^{(k)*}}$ are given. The subproblem is a bilinear quadratic problem due to $\lambda_{nte} d_{nte}$ and $\lambda_{n(i),t,e} g_{ite}$, and the size of the subproblem does not change over the iterations:
\begin{subequations}\label{SP}\begin{align}
    \begin{split}
        &\max_{\Xi^{\mathrm{SP}}} O^{\mathrm{SP}}\!\coloneqq\!\sum_{e\in{\cal E}} \omega_e \Big[
        \hspace{-0.0mm}\sum_{t\in{\cal T}}\sum_{n\in{\cal N}_s}  d_{nte} \pi_{nt}^{(k),*} \hspace{-1mm}
        + \sum_{t\in{\cal T}}\hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}\cup{\cal I}^{\mathrm R}}  \hspace{-4mm}\tau^{\mathrm{e}^{(k),*}}\hspace{-1.5mm} g_{ite}
        \\&\hspace{1.5mm}
        + \hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}} \tau^{\mathrm{c}^{(k),*}} g^{\mathrm{max}}_i
        - \hspace{-0.5mm}\sum_{t\in{\cal T}} \Big( \hspace{-0.5mm}\sum_{n\in{\cal N}_s} \hspace{-1.5mm} \lambda_{nte} d_{nte}\hspace{-1mm}
        -  \hspace{-1.5mm}\sum_{i \in {\cal I}\cup\hat{\cal I}} \hspace{-1.5mm}\lambda_{n(i),t,e} {g}_{ite} \Big)
        \\&\hspace{1.5mm}
        - \sum_{i\in\hat{\cal I}} C^{\mathrm{inv}}_i g^{\mathrm{max}}_{i}
        - \sum_{i\in{\cal I}}\sum_{t\in{\cal T}} C^{\mathrm{g}}_i g_{ite} \Big]
        \label{SP:obj}
    \end{split}
\end{align}
\begin{align}
    & \text{subject to:}\nonumber\\
    & \text{Eqs.}~\eqref{Eq:ramp}\text{\textendash}\eqref{Eq:pdown_ULB}, \eqref{Eq:CC_reform}:\hspace{9.5mm}\text{ML constraints}\\
    & \text{Eqs.}~\eqref{TSO_DSPF}\text{\textendash}\eqref{TSO_LineLim}:\hspace{15mm}\text{LL primal constraints}\label{SP:LLconst}\\
    & \text{Eqs.}~\eqref{DLL_gp}\text{\textendash}\eqref{DLL_theta}:\hspace{15mm}\text{LL dual constraints}\label{SP:DLLconst}\\
    & O^{\mathrm{W}}_{e}=O^{\mathrm{DW}}_{e},~\forall e\in{\cal E}:\hspace{4.7mm}\text{LL strong duality}\label{SP:LLduality}
\end{align}\end{subequations}
where $\Xi^{\mathrm{SP}}\coloneqq \Xi^{\mathrm{U}} \cup \Xi^{\mathrm{W}}$.

\subsubsection{Master problem}
Equation~\eqref{MP:obj} minimizes the objective function of state regulator $s$ over variables of all three levels, where Eqs.~\eqref{MP:ULconstraints}\textendash\eqref{MP:LLconst}  define feasible regions of each level. Eqs.~\eqref{MP:LLconst}\textendash\eqref{MP:LLduality} ensure the LL optimum is achieved. In addition, $k-1$ sets of column-and-cut generation constraints are imposed in Eqs.~\eqref{MP:objcut}\textendash\eqref{MP:strongduality}.
For every C\&CG set, a cut with respect to the ML objective function is generated in Eq.~\eqref{MP:objcut}, with \textit{oracular} knowledge of $\Xi^{\mathrm{U}^{(j)*}}$ given from the subproblem. This allows the master problem to approximate the possible decision-making of the ML and LL problem.
Eqs.~\eqref{MP:DCPF}\textendash\eqref{MP:strongduality} keep auxiliary variables in $\Xi^{\mathrm{Aux}^{(j)}}$ within feasible regions of the LL problem. Note that similarly to the subproblem the master problem is also bilinear quadratic, but the number of auxiliary variables and constraints increase proportionally to the number of iterations:
\begin{subequations}\label{MP}\begin{align}
    \begin{split}
        & \min_{\Xi^{\mathrm{MP}}}  O^{\mathrm{R}}_s = \hspace{-0.5mm} \sum_{t\in{\cal T}}\sum_{e\in{\cal E}} \omega_e \Big[ \sum_{n\in{\cal N}_s} D_{nte} \pi_{nt}
        + \hspace{-2.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}\cup{\cal I}^{\mathrm R}} \hspace{-3mm}g_{ite} \tau^{\mathrm{e}}  \Big]\hspace{-10mm}
        \\&\hspace{20mm}
        + \hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}}  g^{\mathrm{max}}_i  \tau^{\mathrm{c}}
    \end{split}\label{MP:obj}\\
    & \text{subject to:}\nonumber\\
    & \text{Eqs.}~\eqref{Eq:Tariff_budget}\text{\textendash}\eqref{Eq:Radeq}:\hspace{14.8mm}\text{UL constraints}\label{MP:ULconstraints}\\
    & \text{Eqs.}~\eqref{Eq:ramp}\text{\textendash}\eqref{Eq:pdown_ULB}, \eqref{Eq:CC_reform}:\hspace{9.5mm}\text{ML constraints}\label{MP:MLconst}\\
    & \text{Eqs.}~\eqref{TSO_DSPF}\text{\textendash}\eqref{TSO_LineLim}:\hspace{15mm}\text{LL primal constraints}\label{MP:LLconst}
\end{align}
\begin{align}
    & \text{Eqs.}~\eqref{DLL_gp}\text{\textendash}\eqref{DLL_theta}:\hspace{15mm}\text{LL dual constraints}\label{MP:DLLconst}\\
    & O^{\mathrm{W}}_{e}=O^{\mathrm{DW}}_{e},~\forall e\in{\cal E}:\hspace{4.7mm}\text{LL strong duality}\label{MP:LLduality}\\
    & \text{with the C\&CG constraints updated at every iteration $k$ as:}\hspace{15mm}\nonumber\\
    \begin{split}
        & \Bigg[\hspace{-0.5mm}
            \sum_{t\in{\cal T}}\sum_{n\in{\cal N}_s}  d_{nte} \pi_{nt}
            + \sum_{t\in{\cal T}}\hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}\cup{\cal I}^{\mathrm R}}  \tau^{\mathrm{e}} g_{ite}
            \\&\hspace{1mm}
            + \hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}} \tau^{\mathrm{c}} g^{\mathrm{max}}_i
            - \hspace{-0.5mm}\sum_{t\in{\cal T}} \Big( \hspace{-0.5mm}\sum_{n\in{\cal N}_s} \hspace{-1.5mm} \lambda_{nte} d_{nte}\hspace{-1mm}
            -  \hspace{-1.5mm}\sum_{i \in {\cal I}\cup\hat{\cal I}} \hspace{-1.5mm}\lambda_{n(i),t,e} {g}_{ite} \Big)
            \\&\hspace{1mm}
            - \sum_{i\in\hat{\cal I}} C^{\mathrm{inv}}_i g^{\mathrm{max}}_{i}
            - \sum_{i\in{\cal I}}\sum_{t\in{\cal T}} C^{\mathrm{g}}_i g_{ite}
        \hspace{0.5mm}\ge
        \\& \hspace{1.5mm}
            \sum_{t\in{\cal T}}\sum_{n\in{\cal N}_s}  d_{nte} \pi_{nt}
            + \sum_{t\in{\cal T}}\hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}\cup{\cal I}^{\mathrm R}}  \tau^{\mathrm{e}} g_{ite}^{(j)}
            \\&\hspace{1mm}
            + \hspace{-0.5mm}\sum_{i\in\hat{\cal I}^{\mathrm R}} \tau^{\mathrm{c}} g^{\mathrm{max}^{(j-1),*}}_i  \hspace{-2mm}
            - \hspace{-0.5mm}\sum_{t\in{\cal T}} \Big( \hspace{-0.5mm}\sum_{n\in{\cal N}_s} \hspace{-1.5mm} \lambda_{nte}^{(j)} d_{nte}\hspace{-1mm}
            -  \hspace{-1.5mm}\sum_{i \in {\cal I}\cup\hat{\cal I}} \hspace{-1.5mm}\lambda_{n(i),t,e}^{(j)} {g}_{ite}^{(j)} \Big)
            \\&\hspace{1mm}
            - \sum_{i\in\hat{\cal I}} C^{\mathrm{inv}}_i g^{\mathrm{max}^{(j-1),*}}_{i}
            - \sum_{i\in{\cal I}}\sum_{t\in{\cal T}} C^{\mathrm{g}}_i g_{ite}^{(j)}
    \end{split}\label{MP:objcut}\\
    & f_{lte}^{(j)} = \frac{1}{X_l} (\theta_{o(l),te}^{(j)}-\theta_{r(l),te}^{(j)}),~\forall {l \in {\cal L}},~t\in{\cal T},\label{MP:DCPF}\\
    \begin{split}
        &\hspace{-1mm}\sum_{i \in {\cal I}_n\cup\hat{\cal I}_n } \hspace{-1.5mm}{g}_{ite}^{(j)}  + \hspace{-1mm}\sum_{l \vert r(l) = n } \hspace{-1.5mm} f_{lte}^{(j)} \!-\hspace{-1mm}\sum_{l \vert o(l) = n } \hspace{-1.5mm} f_{lte}^{(j)}
        \\&\hspace{20mm}=\begin{dcases}
            &{\hspace{-3mm}d_{nte},\quad\forall n\in{\cal N}_s,~t\in{\cal T},}\\
            &{\hspace{-3mm}D_{nte},\quad\forall n\in{\cal N}\backslash{\cal N}_s,~t\in{\cal T},}
        \end{dcases}
    \end{split}\label{MP:Pbalance}\\
    & 0 \le {g}_{ite}^{(j)} \le \overline{g}^{(j-1),*}_{ite},~\forall i \in {\cal{I}}\cup\hat{\cal{I}},~t\in{\cal T}\!,\label{MP:Tglim}\\
    &-{F}^{\mathrm{max}}_l \le f_{lte}^{(j)} \le {F}^{\mathrm{max}}_l,~\forall {l \in {\cal L}},~t\in{\cal T}\!,\!\!\label{MP:Tllim}\\
    & -\underline{\gamma}_{ite}^{(j)}+\overline{\gamma}_{ite}^{(j)}-\lambda_{n(i),t,e}^{(j)} = - C^{\mathrm g}_i ,\quad\forall i \in {\cal{I}},~t\in{\cal T},\label{MP:DLL_gp}\\
    & \!-\!\underline{\delta}_{lte}^{(j)}\!+\!\overline{\delta}_{lte}^{(j)}\!+\!\xi_{lte}^{(j)} \!-\!\lambda_{r(l),t,e}^{(j)}\!+\!\lambda_{o(l),t,e}^{(j)} \!=\! 0 ,~\forall l \in {\cal{L}},~t\in{\cal T},\label{MP:DLL_fl}\hspace{-1mm}\\
    & -\!\!\!\!\sum_{l \vert o(l) = n}\! \frac{\xi_{lte}^{(j)}}{X_l}~+\!\!\!\!\sum_{l \vert r(l) = n}\! \frac{\xi_{lte}^{(j)}}{X_l} = 0, \quad\forall n \in {\cal{N}},~t\in{\cal T},\label{MP:DLL_theta}\\
    \begin{split}
        & \sum_{t\in{\cal T}} \sum_{i \in {\cal I}\cup\hat{\cal I}} -C_i^{\mathrm{g}} {g}_{ite}^{(j)}
        =  \sum_{t\in {\cal T}} \Big( \sum_{i \in {\cal I}\cup\hat{\cal I}} \! \overline{g}^{(j-1),*}_{ite} \overline{\gamma}_{ite}^{(j)}
        \\&\hspace{10mm}
        - \!\!\sum_{n \in {\cal N}} \!D_{nte} \lambda_{nte}^{(j)}
        + \sum_{l \in {\cal L}} \! \big({F}^{\mathrm{max}}_l \overline{\delta}_{lte}^{(j)}\!+\!{F}^{\mathrm{max}}_l \underline{\delta}_{lte}^{(j)} \big) \Big) \!\!
    \end{split}\label{MP:strongduality}\\
    & \Bigg],\quad\forall e\in{\cal E}, ~j = 2,\dots, k, \nonumber
\end{align}\end{subequations}
where $\Xi^{\mathrm{MP}}\coloneqq (\Xi^{\mathrm{R}}\cup \Xi^{\mathrm{U}} \cup \Xi^{\mathrm{W}})\cup \bigcup_{j=2}^{k}\Xi^{\mathrm{Aux}^{(j)}}$ and $\Xi^{\mathrm{Aux}^{(j)}}\coloneqq \{g_{ite}^{(j)}, \theta_{nte}^{(j)}, f_{lte}^{(j)}, \xi_{lte}^{(j)}, \lambda_{nte}^{(j)}:\text{free};~ \underline{\gamma}_{ite}^{(j)}, \overline{\gamma}_{ite}^{(j)}, \underline{\delta}_{lte}^{(j)}, \overline{\delta}_{lte}^{(j)} \ge 0 \}$.

\section{Case Study}\label{Sec:casestudy}
\begin{figure}[t]
    \centering
    \vspace{-1mm}
        \includegraphics[width=0.82\columnwidth]{Figures/ISO-NE_map_wRPS_updated.pdf}
        \caption{\small An 8-zone ISO-NE  system with RPS goals as of 01/01/2020. \vspace{-3mm}}
        \label{fig:ISONE_8zone}
\end{figure}
\begin{table}[b]
        \centering
        \captionsetup{justification=centering, labelsep=period, font=footnotesize, textfont=sc}
        \caption{Generation capacity by state and type [MW]\vspace{-0mm}}
        \begin{tabular}{L{14mm}  R{7.5mm} R{7.5mm} R{7.5mm} R{7.5mm} R{7.5mm} R{7.5mm}}
            \toprule
            &\multicolumn{1}{c}{ME}&\multicolumn{1}{c}{NH}&\multicolumn{1}{c}{VT}&\multicolumn{1}{c}{MA}&\multicolumn{1}{c}{CT}&\multicolumn{1}{c}{RI}\\
            \midrule
            Wind&221.2&140.5&39.0&681.7&132.5&85.0\\
            Solar&41.4&83.8&306.3&1871.3&464.3&116.7\\
            Nuclear&\textendash&1244.0&620.2&684.7&2116.0&\textendash\\
            Coal&311.8&95.4&\textendash&144.4&744.4&1099.5\\
            Oil&1146.9&400.2&\textendash&1111.7&2212.8&435.0\\
            Natural Gas&3862.7&508.0&\textendash&2249.6&621.4&3491.6\\
            \bottomrule
        \end{tabular}
        \label{table:gencap}
    \vspace{-0mm}
    \vspace{3mm}
        \centering
        \captionsetup{justification=centering, labelsep=period, font=footnotesize, textfont=sc}
        \caption{Investment and operation costs of candidate generators \cite{EIA_cost}\vspace{-0mm}}
        \begin{tabular}{L{29mm}  P{14mm} P{14mm} P{14mm} }
            \toprule
            \multirow{2}{*}{}&\multicolumn{1}{c}{Wind}&\multicolumn{1}{c}{Solar}&\multicolumn{1}{c}{Natural Gas}\\
            \midrule
            Investment Cost, \$/kW&1630&2434&895\\
            Operating Cost, \$/MWh&1.1&0.4&20.0\\
            \bottomrule
        \end{tabular}
        \label{table:genspec}
    \vspace{-0mm}
\end{table}

We use the 8-zone ISO New England (ISO-NE) test system \cite{krishnamurthy20168}, covering six states with the RPS goals as in Fig.~\ref{fig:ISONE_8zone}. Table~\ref{table:gencap} and  Table~\ref{table:genspec} report the current generation mix   and cost characteristics of candidate generators. Expansion decisions of non-strategic states are computed based on the model in \cite{munoz2014engineering}. Similarly to \cite{7464891}, we use a hierarchical clustering algorithm \cite{pitt2000applications} to identify 5 representative days from historical hourly demand and renewable generation forecast data collected from the ISO-NE data library \cite{isone2019final,van2018state}. All fossil-fueled generators are assumed to be committed and the renewable forecast error is zero-mean and has $\sigma_{ite}\!=\!0.1$ for wind and $\sigma_{ite}\!=\!0.2$ for solar, while $\eta\!=\!0.03$. Generation resources in non-strategic  states are assumed to offer their full capacity to the market at their marginal cost. Using the investment recovery period of  10 years and a discount rate of 5\%,  $C^{\mathrm{inv}}_i$ is prorated on a daily basis using a net present value approach. Affine control parameters are set as predefined values $\alpha_{it}=1/\mathrm{card}({\cal I}_s^{\mathrm C}\cup\hat{\cal I}_s^{\mathrm C})$ \cite{jaleeli1992understanding,makarov2009operational}. Demand-side parameters in Eq.~\eqref{eq:flexible_demand} are $N=0.25$ and $M=0.25 D_{nte}+20$ \cite{samadi2010optimal}, while peak TOU hours are set as 13:00-21:00, \cite{Eversource}. All models are implemented using Julia v1.3/JuMP v0.21 and solved by the KNITRO solver v12.1 on an Intel Xeon 2.6 GHz CPU with 124GB of memory. The solution procedure is initialized with 300 different starting points, which improves computational performance and overcomes  potential issues arising from inborn local optimality of the bilinear quadratic problems. The average computation time of the  C\&CG algorithm was 97.4 hours with at most 5 iterations, which is acceptable for long-term planning problems when there is no pressure to instantly enforce the optimal decisions. The input data and code are available in \cite{julia_code_gist}.

\subsection{NH Case: investment, policy, cost decisions}\label{Subsec:NH}
\label{sec:case_nh}

\begin{figure}[t!]
    \centering
    \vspace{-1mm}
        \includegraphics[width=1\columnwidth]{Figures/NH_EPayment_updated.pdf}
        \vspace{-8mm}
        \caption{\small (a) Relationship between the optimal energy-based incentives and electricity tariff and (b) comparison between the consumer payment and the cost of energy-based incentives. \vspace{0mm}}
        \label{fig:NH_EPayment}
    \centering
        \includegraphics[width=1\columnwidth]{Figures/NH_CPayment_updated.pdf}
        \vspace{-8mm}
        \caption{\small (a) Relationship between the optimal capacity-based incentives and electricity tariff and (b) comparison between the consumer payment and the cost of capacity-based incentives. \vspace{-3mm}}
        \label{fig:NH_CPayment}
\end{figure}
\begin{figure}[p]
    \centering
    \vspace{-0mm}
        \includegraphics[width=1.0\columnwidth]{Figures/TOUs.pdf}
        \vspace{-7mm}
        \caption{\small Optimal time-of-use (TOU) electricity tariffs for  (a) energy-based and (b) capacity-based incentives. The darker the color, the higher the value of incentives.  \vspace{1mm}}
        \label{fig:TOUs}
    \centering
    \vspace{-0mm}
        \includegraphics[width=\columnwidth]{Figures/NH_dispatch_joint_1col.pdf}
        \vspace{-7mm}
        \caption{\small Dispatch of generators and demand in NH: (a) before and (b) after implementing RPS goals. \vspace{0mm}}
        \label{fig:NH_dispatch_joint}
    \centering
    \vspace{-1mm}
        \includegraphics[width=\columnwidth]{Figures/ISONE_dispatch_joint_1col.pdf}
        \vspace{-7mm}
        \caption{\small Dispatch of generators and demand in the entire ISO-NE system: (a) before and (b) after implementing RPS goals. \vspace{0mm}}
        \label{fig:ISONE_dispatch_joint}
    \centering
    \vspace{-0mm}
        \includegraphics[width=0.98\columnwidth]{Figures/LMPs_resized.pdf}
        \vspace{-3mm}
        \caption{\small Locational marginal prices (LMPs) in the ISO-NE test system before and after implementing RPS goals. \vspace{-0mm}}
        \label{fig:LMP}
\end{figure}
This section focuses on the NH case because it has the earliest RPS goal among other states. The RPS goal is achieved by installing 602 MW of wind generation resources and requires no enhancement to the fossil-fueled mix (due to a high operating cost). Fig.~\ref{fig:NH_EPayment} and \ref{fig:NH_CPayment} present the optimal incentives and average tariff ($\mathbb{E}_t [\pi_{nt}]$), which is computed for each case based on the hourly TOU values in  Fig.~\ref{fig:TOUs}, to support the RPS goal in NH, as well as itemizes the power utility revenue from both the incentives and tariffs. Both the energy- and capacity-based incentives linearly depend on the average tariff. Indeed, as the average tariff reduces, which in turn leads to a lower consumer payment, greater energy and capacity-based incentives are needed to achieve the same RPS goal. Notably, as Fig.~\ref{fig:NH_EPayment} and \ref{fig:NH_CPayment} show,  the energy-based incentive of  \$20/MWh and the capacity-based incentive of \$600/{kW} yield the same value of the average tariff to achieve the NH RPS goal.

Fig.~\ref{fig:NH_dispatch_joint}\textendash\ref{fig:LMP} compare the average dispatch and LMPs across the entire ISO NE test system before and after achieving RPS goals. Due to the increased wind capacity, NH produces more wind power and  reduces the usage of fossil-fueled resources (e.g. natural gas). Despite this reduction, the amount of power exported from NH to other states increases as shown in Fig.~\ref{fig:NH_dispatch_joint}(b). As a result of these dispatch changes,  LMPs in the system   decrease by $\approx$2.5\% on average as shown in Fig.~\ref{fig:LMP}.

Next, to analyze the sensitivity of the optimal policy decisions in Fig.~\ref{fig:NH_EPayment} and \ref{fig:NH_CPayment}, the current NH RPS goal is compared to the lower and higher targets, which are set to 85\% and 115\% of the original RPS goal. As Fig.~\ref{fig:pi_tau_joint} shows, the linear relationship between the average tariff and the energy- and capacity-based incentives perseveres. Furthermore, the slope in Fig.~\ref{fig:pi_tau_joint} becomes steeper for the high RPS scenario, thus indicating that a higher RPS goal requires either a greater average tariff for consumers or greater incentives.

\subsection{Comparison with the MA and ME cases}\label{Subsec:RPS_Sens}

\begin{figure}[t!]
    \centering
    \includegraphics[scale=1]{Figures/NH_Pi_tau_joint_horizon.pdf}
        \vspace{-7mm}
        \caption{\small Sensitivity of the relationship between the optimal incentives and electricity tariff to an RPS goal. \vspace{2mm}}
        \label{fig:pi_tau_joint}
    \centering
    \includegraphics[scale=1]{Figures/Ret_multi_1col.pdf}
        \vspace{-7mm}
        \caption{\small Sensitivity of the relationship between the optimal incentives and electricity tariff in ME, MA and NH to coal retirements. \vspace{0mm}}
        \label{fig:Ret_multi}
\end{figure}

The results attained in Section~\ref{sec:case_nh} can be extended to the MA and ME cases, which have 41.1 and 100\% RPS goals as in Fig.~\ref{fig:ISONE_8zone}. In the following, we also compared the case with no retirement (base case) and the case with the full retirement of coal-fired resources. As  Fig.~\ref{fig:Ret_multi} shows, the linear relationship between the average tariff and the energy- and capacity-based incentives persevere for all scenarios. However, a more aggressive RPS goal leads to a steeper slope.

The difference between the base and retirement cases manifests itself differently for each state. There is almost no difference in the most aggressive RPS case of ME, where 311.8MW of coal-fired resources can be retired with no implications on the optimal policies. This is due to the fact that the state of ME has a $100\%$ RPS goal, and thus, coal-fired resources would not be dispatched due to a relatively high operating cost once the renewable generation expansion is made. In the case of MA, the retirement will increase the policy incentives needed to  support the RPS goal for the same average tariff.  On the other hand, in the case of NH, coal retirements reduce both the incentives and average tariffs needed to support their modest RPS goal. This is due to the replacement of coal-fired resources in the optimal dispatch with relatively cheap wind and gas producer.

\section{Conclusion}\label{Sec:conclusion}

This paper proposed a TL framework to support a renewable policy decision-making process of a state regulator that internalizes conflicting interests of different stakeholders. The resulting problem is a computationally complex TL bilinear optimization problem and, thus, the C\&CG algorithm is implemented to solve it efficiently. The case study compared the optimal  policies supporting the implementation of RPS goals. The numerical results obtained from the the 8-zone ISO NE test system demonstrated that a more ambitious RPS target requires either an increased electricity tariff (consumer payment) or higher RES incentives (policy support) under the same retirement scenario. Furthermore, our numerical results demonstrate that there is a linear relationship between the RES incentives and  the average  electricity tariff to implement given RPS goals. However, the retirement policy of fossil-fueled generators such as coal-fired power plants can increase or decrease the total cost (incentives and tariffs hikes) of  achieving the same given RPS goals, depending on the current generation mix of the jurisdiction.

\bibliographystyle{IEEEtran}