A Hamiltonian Energy Resource Dispatch Model
A Hamiltonian formulation as an alternative to the Lagrangian formulation of the optimal energy resource dispatch solution when resources are presented with both energy and ramping prices instead of only energy prices.
Several years ago my colleague Sahand Bebhoodi and I wrote a paper [1] on how to optimally dispatch energy resources when they are presented with two prices, one associated with costs as a function of power (quantities in units of Watts and prices in currency units per Joule), and one price for ramping (quantities in Watts per second and prices in currency units per Watt). The motivation for this paper was that dispatching resources using only the conventional energy price (e.g., $/MWh) only works when the marginal units have consistently non-zero marginal costs. If a system has high levels of renewables, the marginal units’ asking prices are theoretically zero a significant fraction of the time. This would seem to defeat the point of having a price signal at all.
Our proposal was to introduce a second price to reveal the cost of ramping. The solution we obtained provided some useful insights, as well as a very simple controller model, both of which are summarized in the next section. However, the solution is not always easy to work with. So in this post, I provide a simpler and potentially more useful model by applying a Legendre transformation to the original solution we found.
Lagrangian Formulation
In the original paper, Sahand and I imagine a situation where there is a secondary price signal available to the distributed resources (in this case a ramping price), and we attempted to answer the question of how one might dispatch resources differently when that second price is considered.
The result was quite interesting and enlightening. We considered a simple situation in which both the power and ramping price are based on quadratic costs functions, and the total cost over a time interval $T$ is simply
\[C(T) = \int_0^T p(q;t) q(t) + r(\dot q;t) \dot q(t) ~ dt\]where $p(q;t) = \frac12 a q(t)$ and $r(\dot q;t) = \frac12 c \dot q(t)$. If we have a constraint on the total energy over the interval, e.g. $\int_0^T q(t)~dt = E$, then we can use a first-order Euler-Lagrange equation
\[\frac{\partial L}{\partial q} - \frac{d}{dt} \frac{\partial L}{\partial \dot q} = 0\]where $L = \frac12 a q^2 + \frac12 c \dot q^2 + \lambda q$ and $\lambda$ is a Lagrange multiplier representing the price of satisfying the energy constraint to find the cost minimal trajectory, which is the solution to the second-order ordinary differential equation
\[c \ddot q - a q = \lambda.\]This equation is easily solved in the complex frequency or $s$-domain, which gives us the solution
\[\hat q(s) = \frac{ Q_0 s^2 + \dot Q_0 s + \frac{\lambda}{c} }{ s(s^2-\omega^2) }\]where $\omega=\sqrt{\frac{a}{c}}$, and $Q_0$ and $\dot Q_0$ are the initial power $q(0)$ and initial ramp rate $\dot q(0)$, respectively.
If we constrain the final power such that $q(T) = Q_T$, then we can solve the following system of equations to obtains the values of $\lambda$ and $\dot Q_0$ such that
\[\begin{bmatrix} \lambda \\ \dot Q_0 \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix}^{-1} \begin{bmatrix} E - D ~ Q_0 \\ Q_T - Q_O ~ \cosh \omega T \end{bmatrix}\]where
- $A = \frac{\sinh\omega T - \omega T}{a\omega}$
- $B = \frac{\cosh\omega T - 1}{\omega^2}$
- $C = \frac{\cosh\omega T - 1}{a}$
- $D = \frac{\sinh\omega T}{\omega}$
Based on this solution, we have the following results. First, when the marginal cost of energy is zero due to high renewables, then the minimum cost dispatch is a quadratric function
\[q(t) = \frac{\lambda}{2c} t^2 + \dot Q_0 t + Q_0\]which is considerably different from the conventional dispatch for energy resources, suggesting that today’s market-based methods of dispatching resources may not be cost-optimal when energy prices are near zero, and ramping prices are non-trivial.
Second, if we do a partial fraction expansion of the $s$-domain solution to the ODE such that $\hat q(s) = \frac{K_1}{s+\omega}+\frac{K_2}{s}+\frac{K_3}{s-\omega}$, then we can divide the dispatch in time domain into three intervals characterized by distinct dominant responses:
- the initial response dominated by the transition $q(t)=K_1 e^{-\omega t}$ from the initial load $Q_0$ to the scheduled load $Q_E = -\frac{\lambda}{a}$;
- the steady response dominated by the scheduled load $q(t)=K_2$; and
- the final response dominated by the transition $q(t)=K_3 e^{\omega t}$ from the scheduled load $Q_E$ to the final load $Q_T$
where $K_1=\frac12 Q_0-\frac{1}{2\omega}\dot Q_0+\frac{\lambda}{2a}$, $K_2=-\frac{\lambda}{a}$, and $K_3=\frac12 Q_0+\frac{1}{2\omega}\dot Q_0+\frac{\lambda}{2a}$.
Hamiltonian Formulation
The formulation of the problem above presents the solution to the optimal resource dispatch problem in phase space, i.e., in terms of the power $q(t)$ and its corresponding ramp rate $\dot q(t)$. This formulation is convenient insofar as it describes the ODE solution quite clearly. However, another formulation exists in terms of the ramping price $r = c \dot q$, which happens to be the conjugate pair variable for the Legendre transformation
\[H(q,r;t) = L(q,\dot q;t) - \dot q \frac{\partial}{\partial \dot q} L(q,\dot q;t)\]because
\[r = \frac{\partial L}{\partial \dot q} = c \dot q.\]This allows us to identify a Hamiltonian function for the system
\[H(q,r;t) = \frac12 a q^2 + \lambda q - \frac1{2c} r^2.\]This function has an important property that emerges from its symplectic structure
\[\left\{ \begin{array}{rll} \dot r &= a q + \lambda &= \frac{\partial H}{\partial q} \\ \dot q &= \frac1c r &= - \frac{\partial H}{\partial r} \end{array} \right.\]i.e., it is constant because
\[\frac{dH}{dt} = \frac{\partial H}{\partial q} \frac{dq}{dt} + \frac{\partial H}{\partial r} \frac{dr}{dt} = \dot r \dot q - \dot q \dot r = 0.\]We now identify two distinct costs (or revenues if negative) associated with these quantities, i.e., the cost of power dispatch and the revenue of ramping dispatch
\[U = \frac12 a q^2 + \lambda q \qquad \mathrm{and} \qquad V = - \frac12 c \dot q^2 = -\frac1{2c} r^2,\]respectively, such that we have
\[L = U - V \qquad \mathrm{and} \qquad H = U + V.\]