\section{ISO budget balancing}
\label{sec:4}

In the previous sections we addressed the problem of stabilizing the wholesale and retail prices to a common value which guarantees the maximization
of social welfare. In the proposed solution, the ISO only decides the next retail price for the consumers according to the energy request during the
previous period. Due to the difference between this prediction and the actual consume, at each time step there exists a gap between the amount of
money to be paid to the producers and the amount of money received from the consumers. This means that the ISO is forced to fill this gap either by
adding money or keeping it, which yields to a non--zero total income or loss, in contrast with the non--for--profit nature of the ISO. The goal of
this section is to provide an effective algorithm to solve this problem, namely, denoting by $I_t, t \geq 0$ the total revenue of the ISO, to steer
it asymptotically to zero.\\
Recalling that $\lambda^r_t$, $\lambda^w_t$ and $x_t$ denote respectively the retail price, the wholesale price and the amount of energy flowing in
the network in the period $[t, t+1)$, we define the interval revenue:
    $$
        i_t \stackrel{\triangle}{=} x_t(\lambda^r_t - \lambda^w_t)
    $$
which is positive if the amount of money received from the consumers exceeds that to be paid to the producer, and negative otherwise. The total
revenue of the ISO is thus simply:
    $$
        I_t \stackrel{\triangle}{=} \sum_{k=1}^ti_k = \sum_{k=1}^tx_k(\lambda^r_k - \lambda^w_k).
    $$
The update law for the retail price proposed in the previous sections was:
    $$
        \lambda^r_{t+1} = \Pi(\lambda^r_t,\, \lambda^w_t) = \lambda^r_t + \gamma(\lambda^w_t - \lambda^r_t)
    $$
which is now modified to:
    $$
        \lambda^r_{t+1} = \Pi(\lambda^r_t,\, \lambda^w_t,\,I_t) = \lambda^r_t + \gamma(\lambda^w_t - \lambda^r_t) - \rho I_t
    $$
where $\rho > 0$ is a real constant whose value is to be decided during the design of the algorithm. Formally, this is equivalent to an integral
control over the difference of the prices weighted according to the energy consumption. Intuitively, what happens is that the larger is the total
income of the ISO, the less will be the retail price in the following periods, since the ISO has a positive amount of money which is to be given back
to the consumers. On the contrary, if the ISO paid in the past to the producers more than what it received by the consumers, then the retail price
will be increased in order to gain more money and fill the cash deficit.\\
The equations of the algorithm in the time period $[t,\,t+1)$ are thus:
    \begin{equation}
        \begin{cases}
            s_t = d_t                                                                                                                 \\
            \lambda^w_{t+1} = \dot{c}(s_t)                                                                                                \\
            \lambda^r_{t+1} = \Pi(\lambda^r_t,\, \lambda^w_t,\,I_t) = \lambda^r_t + \gamma(\lambda^w_{t+1} - \lambda^r_t) - \rho I_t  \\
            d_{t+1} = \dot{v}^{-1}(\lambda^r_{t+1})                                                                                       \\
            I_{t+1} = I_t + i_{t+1} = I_t + d_{t+1}(\lambda^r_{t+1} - \lambda^w_{t+1})                                                           \\
        \end{cases}
    \end{equation}
with initialization $I(0) = 0$.\\
The following result gives a local criterion for the stability of the system.
\begin{theorem}
Assume that, around the equilibrium point for the algorithm $\lambda^*$, it holds:
    \begin{align*}
        &\dot{v}(x) = -ax + b  \\
        &\dot{c}(x) = cx + d.
    \end{align*}
Then the equilibrium point is asymptotically stable if:
    \begin{equation}
        \label{eq:condizione_equiv}
        \begin{cases}
            &   \gamma < \frac{2a}{c+a}                                     \\
            &   \frac{\rho}{2} > \gamma\frac{c+a}{b-d} - \frac{2a}{b-d}
        \end{cases}
    \end{equation}
while a couple of conditions which are only sufficient is:
    \begin{equation}
         \label{eq:condizione_suff}
        \begin{cases}
            &   \gamma < \frac{2a}{c+a} \\
            &   \rho < \gamma\frac{c+a}{b-d}
        \end{cases}.
    \end{equation}
\end{theorem}
\begin{proof}
For simplicity, set $u_t = \lambda^r_t$. By exploiting the assumptions on $\dot{v}$ and $\dot{c}$, we can rewrite the equations of the system as:
    $$
        \begin{cases}
            u_{t+1} = u_t + \gamma(cx_t + d - u_t) - \rho I_t    \\
            I_t = I_{t-1} + (u_t - cx_t - d)x_t                    \\
        \end{cases}
    $$
whence, since $x_t = -\frac{u_t - b}{a}$,
    $$
        \begin{cases}
            u_{t+1} = u_t + \gamma(-\frac{c}{a}(u_t-b) + d - u_t)) - \rho I_t    \\
            I_t = I_{t-1} - (u_t + \frac{c}{a}(u_t - b) - d)\frac{u_t-b}{a}                    \\
        \end{cases}
    .
    $$
The first equation can be rewritten as:
    \begin{align*}
        u_{t+1}
            &   =   u_t - \gamma\frac{c}{a}(u_t - b) + \gamma d - \gamma u_t - \rho I_t \\
            &   =   u_t - \gamma\frac{c}{a}u_t + \gamma \frac{cb}{a} + \gamma d - \gamma u_t - \rho I_t \\
            &   =   u_t(1 - \gamma\frac{c}{a} - \gamma) + \gamma(\frac{bc}{a} + d) - \rho I_t \\
            &   =   f_1(u_t,\,I_t)
        ,
    \end{align*}
while the second equation yields:
    \begin{align*}
        I_{t+1}
            &   =   \left(u_t + \frac{c}{a}u_t - \frac{bc}{a} - d\right)\frac{b - u_t}{a} + I_t   \\
            &   =   \left[u_t\left(1 + \frac{c}{a}\right) - \left(\frac{bc}{a} + d\right)\right]\frac{b - u_t}{a} + I_t  \\
            &   =   -\frac{1}{a}\left(1 + \frac{c}{a}\right)u_t^2 + \frac{b}{a}\left(1 + \frac{c}{a}\right)u_t + \frac{1}{a}\left(\frac{bc}{a} + d\right)u_t - \frac{b}{a}\left(\frac{bc}{a} + d\right) + I_t   \\
            &   =   -\frac{1}{a}\left(1 + \frac{c}{a}\right)u_t^2 + \left(\frac{b}{a}\left(1 + \frac{c}{a}\right) + \frac{1}{a}\left(\frac{bc}{a} + d\right)\right)u_t - \frac{b}{a}\left(\frac{bc}{a} + d\right) + I_t    \\
            &   =   f_2(u_t,\,I_t)
        .
    \end{align*}
The equilibrium point of this system is given, in terms of amount of energy and total revenue of the ISO, by:
    $$
        (x^*,\,I^*) = \left(\frac{b - d}{c + a},\,0\right)
    $$
so that we have:
    $$
        u^* = \lambda^{r*} = -a\frac{b - d}{c+a} + b = \frac{ad + bc}{c + a}
    .
    $$
The next step in the proof is to use the Lyapunov's Linearization Theorem \cite{Fornasini} to prove the claim. In order to do this, we consider the
Jacobian of the system near $(u^*,\,I^*)$, which is given by:
    \begin{align*}
        J(u^*,\,I^*)
        &   =
            \left[
            \begin{array}{cc}
                \frac{\partial f_1}{\partial u}  &  \frac{\partial f_1}{\partial I} \\
                &   \\
                \frac{\partial f_2}{\partial u}  &  \frac{\partial f_2}{\partial I}
            \end{array}
            \right]
           =
            \begin{bmatrix}
                1 - \gamma\frac{c}{a} - \gamma                                                                              &   -\rho   \\
                -\frac{2}{a}(1 + \frac{c}{a})u + \left(\frac{b}{a}(1 + \frac{c}{a}) + \frac{1}{a}(\frac{bc}{a} + d)\right)  &   1
            \end{bmatrix}\vert_{(u^*,\,I^*)}
    .
    \end{align*}
Since:
    \begin{align*}
        &-\frac{2}{a}\left(1 + \frac{c}{a}\right)u + \left(\frac{b}{a}\left(1 + \frac{c}{a}\right) + \frac{1}{a}\left(\frac{bc}{a} + d\right)\right)|_{u^*} =\\
            &   \qquad=   -\frac{2}{a}\frac{c+a}{a}\frac{ad + bc}{c + a}+ \left(\frac{b}{a}\left(1 + \frac{c}{a}\right) + \frac{1}{a}\left(\frac{bc}{a} + d\right)\right)
            \\
            &\qquad
            =   -2\frac{ad + bc}{a^2} + \frac{b(c+a)}{a^2} + \frac{bc + dc}{a^2}
            \\
            &\qquad
            =\frac{b(c+a)}{a^2} - \frac{bc + da}{a^2}
            =
            \frac{1}{a^2}(bc + ba - bc - da)
            \\
            &\qquad
            =
            \frac{b - d}{a}
    \end{align*}
we obtain:
    $$
        J(u^*,\,I^*)
        =
        \begin{bmatrix}
            1 - \frac{\gamma}{a}(c+a)   &   -\rho   \\
            \frac{b-d}{a}               &   1
        \end{bmatrix}
    .
    $$\\
To conclude for the local stability of the algorithm using Lyapunov's Linearization Theorem we must provide conditions so that the two eigenvalues of
$J(u^*,\,I^*)$ are strictly stable, namely inside the unit circle. The characteristic polynomial of the matrix is:
    $$
        P_J(z) = z^2 + z\underbrace{\left(-1 -1 +\frac{\gamma}{a}(c+a)\right)}_{=a_1} + \underbrace{1 - \frac{\gamma}{a}(c+a) -
        \rho\frac{d-b}{a}}_{=a_0}
    ,
    $$
and it is known that in order $P_J(z)$ to have strictly stable roots it is necessary that:
    \begin{align*}
        &   |a_0| < 1   \\
        &   (1 + a_0)^2 - a_1^2 >0
    .
    \end{align*}
Concerning the first condition, it must be true that $a_0 < 1$, so:
    \begin{align*}
        &   1 - \frac{\gamma}{a}(c+a) + \rho\frac{b-d}{a} < 1
        \Rightarrow -\frac{\gamma}{a}(c+a) + \rho\frac{b-d}{a} < 0  \\
        &   \Rightarrow \rho(b-d) < \gamma(c+a)
        \Rightarrow \rho < \gamma\frac{c+a}{b-d}
    \end{align*}
and $a_0 > -1$, which is:
    \begin{align*}
        &   1 - \frac{\gamma}{a}(c+a) + \rho\frac{b-d}{a} > -1
        \Rightarrow \rho\frac{b-d}{a} > -2 + \frac{\gamma}{a}(c+a)\\
        & \Rightarrow \rho > \gamma\frac{c+a}{b-d} - \frac{2a}{b-d}
    \end{align*}
where for both the inequalities we have used $a > 0$ and $b-d>0$, which is needed in order the equilibrium point to be positive, namely physically
feasible. The first condition implies thus:
    \begin{align}
        \label{eq:condizioni_a0}
            &   \rho < \gamma\frac{c+a}{b-d}    \\
            &   \rho > \gamma\frac{c+a}{b-d} - \frac{2a}{b-d}
    \end{align}\\
The second condition implies instead:
    \begin{align*}
        &   \left(2 - \frac{\gamma}{a}(c+a) + \frac{\rho}{a}(b-d)\right)^2 - \left(-2 + \frac{\gamma}{a}(c+a)\right)^2 > 0  \\
        &   \Rightarrow \left(2 - \frac{\gamma}{a}(c+a)\right)^2 + \left(\frac{\rho}{a}(b-d)\right)^2 + \left(2 -
        \frac{\gamma}{a}(c+a)\right)\left(\frac{\rho}{a}(b-d)\right) \\
        &\qquad-\left(-2 + \frac{\gamma}{a}(c+a)\right)^2 > 0 \\
        &   \Rightarrow \frac{\rho^2(b-d)^2}{a^2} + \frac{2}{a^2}\left(2a - \gamma(c+a)\right)\rho(b-d) > 0
    \end{align*}
thus, since, again, $b > d$ and $a> 0$,
    \begin{equation}
        \label{eq:conditioni_a1}
        \frac{\rho}{2} > \gamma\frac{c+a}{b-d} - \frac{2a}{b-d}
    .
    \end{equation}
Notice that this last equation implies, since $\rho > 0$, the second in Equation \eqref{eq:condizioni_a0}, which can be thus discarded. The two remaining
conditions are thus equivalent to the stability of the system. We can moreover notice that if $\gamma < \frac{2a}{c+a}$, then the Right Hand Side of
Equation \eqref{eq:conditioni_a1} is negative, thus the condition is always satisfied. Referring to Fig.~\ref{fig:condizioni}, this means taking $\gamma <
\frac{2a}{c+a}$. A couple of conditions which are only sufficient is thus:
    \begin{equation}
        \begin{cases}
            &   \gamma < \frac{2a}{c+a} \\
            &   \rho < \gamma\frac{c+a}{b-d}
        \end{cases}
        ,
    \end{equation}
and the Theorem is proved.
\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=10cm]{condizioni}
    \end{center}
  \caption{Conditions for the asymptotic stability of the equilibrium point in terms of $\rho$ and $\gamma$. In grey the stability
    region. The point $(\frac{2a}{c+a},\,0)$ is relative to the sufficient condition.}\label{fig:condizioni}
\end{figure}


\subsection{A numerical example}
To give an example of application of the previous theorem, we can consider the following example. We set $c(x) = \frac{1}{2}x^2$ and $v(x) =
20\log(5x)$ for simplicity, which yield $\dot{c}(x) = x$ and $\dot{v}(x) = \frac{20}{x}$. The equilibrium point is $x^* = \sqrt{20}$, for which we
have $a = 1$, $c = \sqrt{20}$, $b = 1$ and $d = 0$. If we set $\gamma = 0.2$ and $\rho = 0.01$ it can be seen that the conditions are satisfied, and
in fact the system is stable and converges to the equilibrium point, as it can be seen in Fig.~\ref{fig:traj_con_cont_int}. In
Fig.~\ref{fig:traj_con_cont_int_zoom} it is depicted the same trajectory zoomed around the equilibrium point, while in Fig.~\ref{fig:ISO_revenue} we
present the behavior of the ISO revenue in time. As it can be seen, at first, due to the small value of $\rho$, the influence of the ISO revenue is
negligible, and the algorithm works in the usual way just trying to steer the retail price and the wholesale price to the same value. Since the
initial energy consumption is on the right with respect to the equilibrium point, the retail price chosen by the ISO is too small, and it has to pay
to the producers more than what it receives from the consumers, making up a shortfall. This grows fast, and despite the fact that the region around
the equilibrium for the prices is reached in a few steps, the ISO maintains the retail prices higher than what it should in order to repay its
deficit. This policy is kept for more than what it should (the integral control is slow), so that the ISO accrues a profit, which is slowly given
back to the consumers. This explains the ``overshoot'' in Fig.~\ref{fig:ISO_revenue} and the fact that, in Fig.~\ref{fig:traj_con_cont_int_zoom}, we
can notice that the price is at first higher than the equilibrium and than slightly lower.

\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=10cm]{ISO_revenue}
    \end{center}
  \caption{Trajectory of the total income/shortfall of the ISO.}\label{fig:ISO_revenue}
\end{figure}

\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=10cm]{traj_con_cont_int}
    \end{center}
  \caption{Trajectory of the energy request vs the retail price.}\label{fig:traj_con_cont_int}
\end{figure}

\begin{figure}[ht!]
    \begin{center}
        \includegraphics[width=10cm]{traj_con_cont_int_zoom}
    \end{center}
  \caption{Trajectory of the energy request vs the retail price. Zoom around the equilibrium point.}\label{fig:traj_con_cont_int_zoom}
\end{figure}
\end{proof}
