\section{Stability Analysis}
\label{sec:3}

The goal of this section is to present a simple stabilizing algorithm, based on the main outcome proposed in \cite{Sanjoy2}, which leads to the same
value both the retail and wholesale prices, value which is exactly the optimal value in the sense of the global benefit. The steady assumption is
that the consumers' behavior is static, namely the value function is time--independent. The proof is very simple and is based on the linearization
method for discrete time systems of the Lyapunov Stability Analysis.
\begin{theorem}
Suppose that there exist a maximum value of the demand
$d_{max} > 0$, and for each trajectory satisfying the model equations \ref{modello} consider the demand $d_{t} \in [0, d_{max}]$.
Let the pricing function $\Pi$ be defined as:
\begin{equation}
\Pi(\lambda^{r},\lambda^{w}) = \lambda^{r} + \gamma (\lambda^{w}-\lambda^{r}) = \gamma \lambda^{w} + (1-\gamma)\lambda^{r}
\end{equation}
where $\gamma > 0$. Then for sufficiently small $\gamma$, the pricing function $\Pi$ stabilizes the price's dynamic:
\begin{equation}
\lambda^{r}_{t+1} = \Pi(\lambda^{r}_{t},\lambda^{w}_{t}) = \lambda^{r}_{t} + \gamma (\lambda^{w}_{t}-\lambda^{r}_{t})
\end{equation}
in the sense that $\lambda^{r}_{t}$ and $\lambda^{w}_{t}$ converge to the same value $\lambda^{r*} = \lambda^{w*} = \lambda^{*}$,
and the corresponding supply and demand converge also to equilibrium: $s^{*} = d^{*} = x^{*}$.\\
Let's assume that the equations for the dynamic of the retail and wholesale prices are given by:
\begin{eqnarray}
\lambda^{w} = \dot c (x) = c \cdot x, \label{cpunto} \\
\lambda^{r} = \dot v (x) = a(x), \label{vpunto}
\end{eqnarray}
where $c>0$, because the cost function is strictly convex, while $a(x) > 0$ and its derivative $\dot a(x) < 0$, because the value function is strictly concave and so its
derivative is monotonically decreasing. In such a context the interval of admissible values for the equilibrium's stability
is $\gamma \in [0, \gamma^{*})$, where:
\begin{equation}
\gamma^{*} = \frac{2 \, \dot a(x^{*})}{\dot a(x^{*}) - c} \label{gamma}.
\end{equation}
\end{theorem}
\begin{proof}
Let's start giving an expression of the dynamics of the retail price as function of the same variable $\lambda^{r}_{t}$ only.
In particular we have:
\begin{eqnarray}
\lambda^{r}_{t+1} &=& \lambda^{r}_{t} + \gamma (\lambda^{w}_{t}-\lambda^{r}_{t}) \nonumber \\
            &=& \lambda^{r}_{t}(1-\gamma) + \gamma \lambda^{w}_{t} \nonumber \\
                          &=& \lambda^{r}_{t}(1-\gamma) + \gamma \, c \, x(t) \nonumber \\
                          &=& \lambda^{r}_{t}(1-\gamma) + \gamma \, c \,a^{-1}(\lambda^{r}_{t}) \nonumber \\
                          &=& f(\lambda^{r}_{t}) \label{f}.
\end{eqnarray}
We note that $f(\cdot)$ is not a linear function of the variable $\lambda^{r}_{t}$, which is continuous with its derivatives, and $\lambda^{*}$ is the point of equilibrium.
We calculate now the derivative of $f(\cdot)$ w.r.t.\  $\lambda^{r}$, (the Jacobian matrix in a multivariable case), writing for simplicity $\lambda$ in place of
$\lambda^{r}$:
\begin{eqnarray}
\left.\frac{\partial f}{\partial \lambda}\right|_{\lambda^{*}} &=& \left.1 - \gamma + \gamma \, c \, \frac{1}{\dot a (a^{-1}(\lambda))} \right|_{\lambda^{*}} \nonumber \\
                                                                                         &=& 1 - \gamma + \gamma \, c \, \frac{1}{\dot a (x^{*)}} \label{df}.
\end{eqnarray}
To apply the linearization method for discrete time systems of Lyapunov Analysis (\cite{Fornasini} Chap.4 pg. 179),
we now have to check if the modulus of the eigenvalues of the Jacobian Matrix, which is in this case given only by the value in \eqref{df},
is smaller than 1. So we impose:
\[
|1 - \gamma + \gamma \, c \, \frac{1}{\dot a (x^{*)}}|<1
\]
which is equivalent to:
\begin{eqnarray}
-1 &<& 1 - \gamma + \gamma \, c \, \frac{1}{\dot a (x^{*)}} < 1\nonumber \\
-2 &<& - \gamma + \gamma \, c \, \frac{1}{\dot a (x^{*)}} < 0 \nonumber \\
-2 &<& \gamma \, ( -1 + c \, \frac{1}{\dot a (x^{*)}}) < 0 \label{mod}.
\end{eqnarray}
The second bound in \eqref{mod} is always verified, because the term within brackets is negative by assumption.\\
On the other side, considering the first bound in \eqref{mod}, we find:
\begin{equation}
\gamma < \frac{2 \, \dot a(x^{*})}{\dot a(x^{*}) - c}
\end{equation}
which is exactly the equation \eqref{gamma} of the Theorem.
\qed
\end{proof}
The dynamic of prices, subject to the presented stabilizing algorithm, can be represented
with the scheme in Figure \ref{fig:alg_stab}.
\begin{figure} [!hh]
\centering
\epsfig{file=algoritmo_base_schema.pdf, width=0.65\columnwidth}
\caption{Block diagram describing the stabilizing pricing algorithm. \label{fig:alg_stab}}
\end{figure}

\begin{remark}
We have found an interesting upper bound for the control parameter $\gamma$ in the dynamics of prices that ensures the stability of the system. In
the limiting case of choosing $\gamma=\gamma^{*}$, which corresponds to an eigenvalue on the unit circle, one should check the value of the higher
order derivatives of the function $f(\cdot)$ of equation \eqref{f}, to have some other information about the stability of the equilibrium. There are
in fact some criteria for first order discrete time systems, that allow to verify the asymptotic stability of the equilibrium in the case of an
eigenvalue with modulus equal to 1 (\cite{Fornasini} chap.4 pg. 181). A special remark has to be made regarding the application of this analytic
result. In fact, as already mentioned, in real applications we do not suppose to know the behaviour of the costumers, and so the value function is
not available. So the exact value of the admissible bound for $\gamma$ is not computable, but anyway the theory demonstrates that an exact bound
exists. The result could anyway taken into account if at least a rough approximation of the value function is known. Clearly, the rougher the
information of the value function, the more conservative the bound on $\gamma$ will be.
\end{remark}\\
The dynamics of the system has been simulated in \textsc{Matlab},
assuming that the cost and value functions of producers and consumers are:
\begin{eqnarray}
c(x) &=& 0.5 \, x^{2}, \\
v(x) &=& 20 \, \log(5x +1),
\end{eqnarray}
whose derivatives are:
\begin{eqnarray*}
\dot c(x) &=& x, \\
\dot v(x) &=& 100 \,(5x +1)^{-1}.
\end{eqnarray*}
Therefore the functions respect the hypothesis of the theorem and we can calculate analytically the superior bound for the
control parameter $\gamma$, which is equal to $\gamma^{*} = 0.9776$.\\
First of all we simulated the dynamics of the system with two
different values of the parameter, first with $\gamma = 0.1$ and then with $\gamma = 0.8$.
The stable dynamics of the retail prices are represented in Figures \ref{fig:gamma_fisso_01} and
\ref{fig:gamma_fisso_08}. 
The increasing line is the function $\dot c(x)$, which gives the wholesale prices
w.r.t.\ the energy units consumption.
Analogously, the decreasing curve is the $\dot v(x)$, which gives the retail prices.
The red line represents the trajectory of the retail prices,
which are, at every step, a convex combination of the wholesale and 
retail prices at the previous time step.\\
In Figures \ref{fig:diff_prezzi_01} and \ref{fig:diff_prezzi_08} the trend of the difference between the retail and wholesale prices, in both cases, is represented.
We notice that, for a small value of the parameter,
$\gamma$ the price difference smoothly converges to zero, while, for a value which is next to the superior bound of the parameter $\gamma^{*}$,
the trend begins to oscillate, before settling to the equilibrium.

\begin{figure} [!h]
\begin{center}
\epsfig{file=traj_gamma_fisso_01_x08.pdf, width=0.65\columnwidth}
\end{center}
\caption{Dynamics of prices with $\gamma = 0.1$ \label{fig:gamma_fisso_01}}
\end{figure}

\begin{figure} [!h]
\begin{center}
\epsfig{file=traj_gamma_fisso_08_x08.pdf, width=0.7\columnwidth}
\end{center}
\caption{Dynamics of prices with $\gamma = 0.8$ \label{fig:gamma_fisso_08}}
\end{figure}

\begin{figure} [!h]
\begin{center}
\epsfig{file=differenza_lambdarelambdaw_gamma_fisso_01_x08.pdf, width=0.7\columnwidth}
\end{center}
\caption{Discrepancies between retail and wholesale prices, with $\gamma = 0.1$ \label{fig:diff_prezzi_01}}
\end{figure}

\begin{figure} [!h]
\begin{center}
\epsfig{file=differenza_lambdarelambdaw_gamma_fisso_08_x08.pdf, width=0.7\columnwidth}
\end{center}
\caption{Discrepancies between retail and wholesale prices, with $\gamma = 0.8$ \label{fig:diff_prezzi_08}}
\end{figure}
In Figures \ref{fig:stab} and \ref{fig:instab} it is depicted respectivly the dynamics of the prices when the parameter $\gamma$ is chosen slightly smaller 
or bigger than the critical value. As it is expected, in the first case
the system is stable, but exhibits big oscillations, while in the second
it exits from the linear region, namely the equilibrium point is unstable.
\begin{figure} [!h]
\begin{center}
\epsfig{file=traj_gamma_limite_stabile_x010.pdf, width=0.7\columnwidth}
\end{center}
\caption{Dynamics of prices with $\gamma = 0.98 \,\gamma^{*} $ \label{fig:stab}}
\end{figure}

\begin{figure} [!h]
\begin{center}
\epsfig{file=traj_gamma_limite_instabile_x010.pdf, width=0.7\columnwidth}
\end{center}
\caption{Dynamics of prices with $\gamma = 1.02 \,\gamma^{*} $ \label{fig:instab}}
\end{figure}

\subsection{Estimation of the value function and convergence rate maximization}

In the previous section we provided a control algorithm which forces the retail price to converge to the equilibrium, i.e.\ the value corresponding
to the consumption of $x^*$ units of energy, such that $\dot{v}(x^{*})=\dot{c}(x^{*})$. The goal of this section is to provide a design algorithm to
maximize the rate of convergence by a clever choice of $\gamma$. Recalling that the derivative of the cost function for the suppliers $\dot{c}(x)$ is
supposed to be known, it is trivial that, under the usual hypothesis on the concavity of the value function, there exists an optimal value of
$\gamma$.
Moreover, if we knew the  exact expression of the $\dot{v}(x)$ function, we would also be able to compute the intersection of the functions and therefore to reach the equilibrium point in one step. In other words, there always exists a $\bar{\gamma}$ such that the ISO controls the system to the optimal price in a single time step using the proposed algorithm.\\
In fact:
\begin{eqnarray*}
\lambda_{r}(t+1) & = & \lambda_{r}(t)+\gamma\left(\lambda_{w}(t)-\lambda_{r}(t)\right)\\
\\\exists\bar{\gamma} & : & \bar{\gamma}=\frac{\lambda_{r,ott}-\lambda_{r}(t)}{\lambda_{w}(t)-\lambda_{r}(t)};
\end{eqnarray*}
and $0<\bar{\gamma}<1$.\\
Unfortunately, in most practical applications the function $\dot{v}(x)$ is unknown, so it is not possible to find the optimal $\bar{\gamma}$ without an a priori knowledge. Still, given partial information about the customers' response to the previously proposed prices on earlier steps, it is possible to get an estimate of the real curve by interpolation with a parameterized curve.\\
At the first step, only one point of the $\dot{v}$ function is known. In this case there is no other way to proceed  than running the algorithm with a sufficient small value of $\gamma$. At the second step, it is possible to interpolate the two points with a line - which is unique - and take $\dot{c}(\bar{x})$ as the next retail price, where $\bar{x}$ is such that $\dot{v}(\bar{x})=\dot{c}(\bar{x})$. At this stage, three points of the value function derivative can be used, and therefore sthere are several possible fitting curves, such as: 1)  a line, 2) a parabola and 3) an hyperbola, all achieved with Ordinary Least Squares (OLS).
\begin{enumerate}

\item \textit{Line}: as a first attempt, we approximated the $\dot{v}$ function with a decreasing line ($f(x)=-ax+b$ with $a>0,\,b>0$). The resulting fit is quite different from the original curve, except for a limited region around the equilibrium point.
Moreover, if all of the previous data are used at every step, it takes a very long time to converge to the equilibrium. The effect of taking into account all the previous points without considering a suitable forgetting factor is to slow down the algorithm, as it is equivalent to weighting too much the initial data.
Therefore it is expedient to fit the latest points only. In this case the convergence rate to the equilibrium point is much higher.

\item \textit{Parabola}: there is a basic problem when the $\dot{v}$ function is approximated with a convex parabola, namely $f(x)=ax^2+x+c$ with $a>0$. The curve, for $x \rightarrow \infty$, diverges, while $\dot{v}(\infty)=0$. It is easy to notice that this fact can lead to two different scenarios. If there are two different intersections between $\dot{c}$ and $\dot{v}$, it is sufficient to take that corresponding to the smaller amount of energy, which is the left one. Otherwise, there can be no intersection at all, in which case the algorithm stops.

\item \textit{Hyperbola}: the approximation throught an equilateral hyperbola ($f(x)=a/x$ with $a>0$) is the most performant among the presented functions. The hyperbola shares many characteristics with the original $\dot{v}$ function, as it is strictly positive, monotonically decresing, convex, with $f(0)\gg 0$ and $f(\infty)=0$.
\end{enumerate}

\begin{remark}
As in this case all the available points belong exactly to the original $\dot{v}$ function, a spline interpolation could have been used. Nevertheless, in most pratical applications the data are affected by some kind of noise, and because of this the spline interpolation could lead to a completely different function.
\end{remark}

\subsection{Numerical results}
We now present the simulation results with the proposed interpolation functions. We set $c(x)=0.5\,x^2 $ as the cost function, and $v(x)=40 \,\log(5x+1)$ as
the value function. This yields to $\dot{c}(x)=x$ and $\dot{v}(x)=200 \,(5x+1)^{-1}$ respectively. At the first step we do not have any a priori information or
bound for $\gamma$, therefore the best we can do is to choose a sufficiently small value for the parameter (in this case, $\gamma=0.1$).\\
\begin{figure} [!h]
\begin{center}
\epsfig{file=interpolazione_con_retta.pdf, width=0.65\columnwidth}
\end{center}
\caption{Line interpolation\label{fig:line}}
\end{figure}
\begin{figure} [!h]
\begin{center}
\epsfig{file=interpolazione_con_retta_ultimi.pdf, width=0.65\columnwidth}
\end{center}
\caption{Line interpolation (last points only)\label{fig:line_ult}}
\end{figure}
\begin{figure} [!h]
\begin{center}
\epsfig{file=interpolazione_con_parabola.pdf, width=0.65\columnwidth}
\end{center}
\caption{Parabola interpolation \label{fig:parabola}}
\end{figure}
\begin{figure} [!h]
\begin{center}
\epsfig{file=interpolazione_con_iperbole.pdf, width=0.65\columnwidth}
\end{center}
\caption{Hyperbola interpolation\label{fig:iperbole}}
\end{figure}\\
The starting price for the simulations is the one that let the customers consume exactly $x(0)=10$ units of energy. As we can see in Figure
\ref{fig:line}, the presence of the first point prevents the line from converging to the equilibrium. The convergence rate is much higher when only
the last four points are used for interpolation (as shown in Figure \ref{fig:line_ult}, as the algorithm is allowed not to take into permanent
account the initial values. The parabola interpolation (Figure \ref{fig:parabola}) is quite fast, but it is not possible to guarantee the convergence
for every initial value. As we can see from Figure \ref{fig:iperbole}, the hyperbolic interpolation is the most accurate, as it provides both high
performance and convergence guarantee.\\
In Table \ref{tab1} we show how many steps are needed to reach a retail price within a range of $5\%$ (or $2\%$) to the optimal retail price.
The starting point is the price that let the customers consume exactly $x(0)=10$ units of energy.\\
\begin{table} [!h]
\begin{centering}
\begin{tabular}{|c|c|c|}
\hline
 & $5\%$ & $2\%$\tabularnewline
\hline
\hline
Line & $3$ & $4$\tabularnewline
\hline
Line (only with the last points) & $3$ & $4$\tabularnewline
\hline
Parabola & $3$ & $3$\tabularnewline
\hline
Hyperbola & $3$ & $3$\tabularnewline
\hline
No interpolation (fixed $\gamma$) & $9$ & $13$\tabularnewline
\hline
\end{tabular}\caption{Number of steps to reach $5\%$ and $2\%$ to the equilibrium ($x(0)=10$).}
\label{tab1}
\par\end{centering}
\end{table}\\
In Table \ref{tab2} are shown the results obtained from a starting price that let the customers consume exactly $x(0)=2$ units of energy.
\begin{table} [!h]
\begin{centering}
\begin{tabular}{|c|c|c|}
\hline
 & $5\%$ & $2\%$\tabularnewline
\hline
\hline
Line & $4$ & $7$\tabularnewline
\hline
Line (only with the last points)  & $4$ & $6$\tabularnewline
\hline
Parabola & $\infty$ & $\infty$ \tabularnewline
\hline
Hyperbola & $3$ & $4$\tabularnewline
\hline
No interpolation (fixed $\gamma$) & $16$ & $20$\tabularnewline
\hline
\end{tabular}\caption{Number of steps to reach $5\%$ and $2\%$ to the equilibrium ($x(0)=2$).}
\label{tab2}
\par\end{centering}
\end{table}\\
Similar results can be obtained using a different value function, for example $v(x)=40 \log(5x+1)$, as shown in Table  \ref{tab3} (where the initial energy consumption is guaranteed to be $x(0)=10$ units) and Table \ref{tab4} (initial energy consumption set to $x(0)=2$ units).
\begin{table}  [!h]
\begin{centering}
\begin{tabular}{|c|c|c|}
\hline
 & $5\%$ & $2\%$\tabularnewline
\hline
\hline
Line & $3$ & $34$\tabularnewline
\hline
Line (only with the last points)& $3$ & $3$\tabularnewline
\hline
Parabola & $3$ & $3$\tabularnewline
\hline
Hyperbola & $3$ & $3$\tabularnewline
\hline
No interpolation (fixed $\gamma$) & $8$ & $12$\tabularnewline
\hline
\end{tabular}\caption{Number of steps to reach $5\%$ and $2\%$ to the equilibrium ($x(0)=10$).}
\label{tab3}
\par\end{centering}
\end{table}
\begin{table}  [!h]
\begin{centering}
\begin{tabular}{|c|c|c|}
\hline
 & $5\%$ & $2\%$\tabularnewline
\hline
\hline
Line & $10$ & $34$\tabularnewline
\hline
Line (only with the last points) & $6$ & $7$\tabularnewline
\hline
Parabola & $\infty$ & $\infty$ \tabularnewline
\hline
Hyperbola & $3$ & $10$\tabularnewline
\hline
No interpolation (fixed $\gamma$) & $20$ & $24$\tabularnewline
\hline
\end{tabular}\caption{Number of steps to reach $5\%$ and $2\%$ to the equilibrium ($x(0)=2$).}
\label{tab4}
\par\end{centering}
\end{table}
