Adaptive Real-Time Grid Operation via Online Feedback Optimization with Sensitivity Estimation

For each bus $i$ of a $n$-bus power system we define the voltage magnitude as $v_{i}\in\mathbb{R}$, the active and reactive power as $p_{i}\in\mathbb{R}$ and $q_{i}\in\mathbb{R}$, respectively. We obtain the vectors $v$, $p$, and $q$ of dimension $n$ by stacking the individual bus quantities, i.e., $v=[v_{1},\dots,v_{n}]^{T}$. We define the control input vector $u\in\mathbb{R}^{n_{u}}$ consisting of all the controllable resources (e.g. active and reactive generation and flexible loads in $p$ and $q$, slack bus voltage magnitude $v_{1}$ through tap changers); the output vector $y$ (e.g. voltage magnitude elements in $v$) with all the quantities that we measure and want to control through the inputs; and the disturbance vector $d$ with all uncontrollable power injections (e.g. conventional consumption loads in $p$ and $q$). The grid admittance matrix and the power flow equations allow to define an input-output map $\mathdutchcal{h}(\cdot)$ that characterizes the output $y$ as a non-linear function of $u$ and $d$:

The input-output map $\mathdutchcal{h}(\cdot)$ is not typically available in closed form, since in general it is not possible to derive an analytical expression of $v$ (in $y$) as a function of $p$ and $q$ (in $u$ and $d$) using the power flow equations [1]. Yet, the local existence of a continuous differentiable map $\mathdutchcal{h}(\cdot)$ can be guaranteed by the implicit function theorem [20].

The operation of a power grid consists of deciding the input $u_{t}$ at each time instant $t$. An AC-OPF allows to formulate this decision process as an optimization problem:

where $f(u)$ is the operational cost on the input $u$; $g(y)$ is a penalty function to enforce some grid specification on the output $y$, e.g., voltage limits; $\mathcal{U}_{t}$ is the time-varying set of admissible inputs that defines the operational constraints on $u_{t}$, e.g., power limits $\mathcal{U}_{t}=\{u\,|\,\underline{u}_{t}<u<\overline{u}_{t}\}$; and $d_{t}$ is the disturbance value at time $t$, e.g., uncontrollable loads or non-dispatchable generation. The nonlinear input-output model (1) in (2) relates the outputs to the chosen input.

Optimal real-time decision making consists of first taking measurements $d_{t}$; then, solving the AC-OPF problem (2), and finally applying the solution $u^{*}_{t}$ to the system. Then, this is repeated at the next time step $t+1$.

Solving AC-OPF problems (2) to determine the set-points of power resources is a compelling and valuable tool for grid operators, but it comes with some drawbacks: First, the full nonlinear model of the grid $\mathdutchcal{h}(u,d)$ is needed. Second, solving the AC-OPF (2) can be computationally expensive, which may jeopardize its use for real-time grid operation. This can be circumvented by linearizing the map $\mathdutchcal{h}(\cdot)$ in (1) at an operating point [21, 22, 1], e.g., the zero-injection point $(u_{\text{op}},d_{\text{op}})=(0,0)$, to obtain the approximation

where $y_{0}$ is an offset representing the output value when $u=d=0$, e.g., $1$ p.u. for all voltage magnitudes. The matrices $H_{0}=\nabla_{u}\mathdutchcal{h}(u,d)|_{(u_{\text{op}},d_{\text{op}})}$ and $D_{0}=\nabla_{d}\mathdutchcal{h}(u,d)|_{(u_{\text{op}},d_{\text{op}})}$ are evaluated at the operating point, and represent the sensitivities of the output with respect to changes in the input $u$ and disturbance $d$, respectively. This linear approximation (3) can substitute the nonlinear map $\mathdutchcal{h}(\cdot)$ in the AC-OPF (2) to get

Solving the AC-OPF with linear power flow approximation (4) is computationally efficient and could be employed in real-time operation. However, this approach does not take advantage of output measurements $y_{t}$, since it only feeds $d_{t}$ through the inaccurate linear model (3). Hence, such a feedforward approach introduces a model-mismatch that can cause a performance degradation, and even lead to constraint violations, e.g., under and overvoltages.

Instead, OFO is a novel approach [5, 6, 4] that uses $y_{t}$ as feedback to achieve a safer grid operation and track the solution of the AC-OPF (2) under time-varying conditions. For that, OFO turns a standard optimization algorithm, in our case projected gradient decent [23], into a feedback controller that takes the grid output measurements $y_{t}$, instead of computing the output $y_{t}$ via the grid model (1) or the linearized one (3). Projected gradient decent consists of a gradient step and a projection: First, we compute the gradient of the cost function in (4):

To minimize the operational cost, the current input $u_{t}$ is pushed along the direction of the negative gradient with a step size $\alpha$, and then it is projected onto the feasible space $\mathcal{U}_{t}$ to enforce the operational constraints on the input, i.e.,

where $\Pi_{\mathcal{U}}\big{[}u]=\arg\min_{z\in\mathcal{U}}\lVert u-z\rVert_{2}^{2}$ is the projection of $u$ onto $\mathcal{U}$, which is typically easy to evaluate for power grid operation [6], especially if $\mathcal{U}_{t}=\{u\,|\,\underline{u}_{t}\leq u\leq\overline{u}_{t}\}$ is a box constraint.

Due to the non-linearity of $\mathdutchcal{h}(u,d)$, the true sensitivity $\nabla_{u}\mathdutchcal{h}(u,d)$ depends on the values of $u$ and $d$. The temporal variation of the disturbance $d_{t}$ and the input $u_{t}$, e.g., due to applying the OFO controller (5) in the input case, produces a time-varying sensitivity $H_{t}=\nabla_{u}\mathdutchcal{h}(u_{t},d_{t})$. Instead of learning the dependency on $u$ and $d$, we model a time-varying sensitivity $H_{t}$ with the following random process:

where $h=\text{vec}(H)$ is the column-wise vector representation of the sensitivity matrix $H$, $\Delta u_{t-1}=u_{t}-u_{t-1}$ denotes a change of the input $u$, and $\omega_{p,t}\sim\mathcal{N}(0,\Sigma_{p,t})$ is a Gaussian process noise with covariance $\Sigma_{p,t}=\Sigma_{p_{1}}+\Sigma_{p_{2}}\lVert\Delta u_{t}\rVert_{2}^{2}$, that represents how the sensitivity changes over time. We make the part $\Sigma_{p_{2}}$ of the process noise proportional to $\lVert\Delta u_{t}\rVert_{2}$, since a large $\Delta u_{t}$ can trigger a larger change in the true sensitivity $\nabla_{u}\mathdutchcal{h}(u,d)$ that depends on $u$, and the part $\Sigma_{p_{1}}$ independent of $\Delta u_{t}$ to account for a uncontrolled random change $\Delta d_{t}=d_{t+1}-d_{t}$ that can affect the sensitivity as well.

Next, to derive a measurement equation for the sensitivity $H_{t}$, consider the first-order Taylor approximation of $y_{t}$

At each time $t$, we measure $y_{t}$, and compute the variation $\Delta y_{t-1}=y_{t}-y_{t-1}$. Based on the Taylor approximation (7), we treat this variation $\Delta y_{t-1}$ as a noisy linear measurement of $H_{t-1}$ through a measurement model that depends on $\Delta u_{t-1}$:

where $U_{\Delta,t}=\Delta u_{t}^{T}\otimes\mathbbm{1}$, with the Kronecker product $\otimes$, and $\omega_{m,t}\sim\mathcal{N}(0,\Sigma_{m,t})$ is a Gaussian measurement noise with covariance $\Sigma_{m,t}=\Sigma_{m_{1}}+\Sigma_{m_{2}}\lVert\Delta u_{t}\rVert_{2}^{2}+\Sigma_{m_{3}}\lVert\Delta u_{t}\rVert_{2}^{4}$. Again, the part $\Sigma_{m_{1}}$ independent of $\Delta u_{t}$ in the measurement noise represents the effect of an uncontrolled random disturbance change $\Delta d_{t}$, while the other parts $\Sigma_{m_{2}}$ and $\Sigma_{m_{3}}$ encapsulate the second-order error of the Taylor approximation (7).

To update the sensitivity estimate $\hat{h}_{t}$, we combine the information given by the previous sensitivity estimate $\hat{h}_{t-1}=\text{vec}(\hat{H}_{t-1})$, and the measurements $\Delta y_{t-1}$ (8). We compute the new sensitivity estimate $\hat{h}_{t}$ through a Bayesian update represented in the following least-squares problem [17, 18]:

where $\Sigma_{t}$ is the covariance matrix representing the uncertainty of the sensitivity estimate $\hat{h}_{t}$, and $\lVert x\rVert_{A}^{2}=x^{T}Ax$ is the norm of $x$ with respect to a positive definite matrix $A$. The resulting recursive estimation can be expressed as a Kalman filter [24]:

where $\mathbbm{1}$ is the identity matrix, and $K_{t}=\Sigma_{t}U_{\Delta,t}^{T}(\Sigma_{m,t}+U_{\Delta,t}\Sigma_{t}U_{\Delta,t}^{T})^{-1}$ is the Kalman gain, which is well defined for an invertible $\Sigma_{m,t}$, see later Assumption 1.

To learn the time-varying sensitivity $H_{t}$, we need to capture enough information via the measurement equation (8), i.e, we need to use different $\Delta u$ to explore different reactions $\Delta y$ and infer different elements of $H_{t}$ from them. This can be formalized via the persistency of excitation condition [25]: $\Delta u_{t}$ is persistently exciting if there exists a time span $T>0$, such that for all $t>0$, the matrix formed by columns $\Delta u_{t+i}$ for $i\in\{0,\dots,T\}$ has full rank, i.e., $\text{rank}(\Delta u_{t},\dots,\Delta u_{t+T})=n_{u}$. To achieve persistency of excitation, we perturb the OFO step (5) with $\omega_{u,t}\in\mathbb{R}^{n_{u}}$, a bounded zero-mean white noise with independent and identically distributed elements with standard deviation $\sigma_{u}$, e.g., a truncated Gaussian distribution. As a result, we obtain the following persistently exciting OFO with estimated sensitivity $\hat{H}_{t}$:

The resulting interconnected OFO, sensitivity learning and power grid is represented in the block diagram in Figure 1. At each time $t$, a complete loop of the online optimization with sensitivity estimation can be represented as:

In this section we analyze the convergence of the estimated sensitivity $\hat{H}_{t}$ produced by the sensitivity learning (9), and the input $u_{t}$ produced by the OFO (10), towards the true sensitivity $H_{t}$ and the solution $u^{*}_{t}$ of the AC-OPF (2), respectively. We certify this convergence assuming that the true sensitivity $H_{t}$ behaves according to the simplified dynamic process (6) and satisfies the linear measurements equation (8); and that the projected gradient descent used in (10) is a strongly monotone and Lipschitz continuous operator:

The continuous differentiability of $f(\cdot)$ and $g(\cdot)$ is common for typical cost functions in power systems, e.g., linear or quadratic $f(\cdot)$, and quadratic penalty functions like $g(\cdot)=\max(0,\cdot)^{2}$. For strongly convex and Lipschitz smooth cost functions $f(\cdot)$, the strong monotonicity and Lipschitz continuity of the gradient operator $F_{t}(\cdot)$ holds in certain regions around nominal operating points [10]. In particular, it would hold if using a usual linear approximation for the input-output map (1) [8]. Since $u$ and $d$ are restricted by the grid physical limits, e.g., power ratings, the upper bound of $\lVert\Delta u_{t}\rVert_{2}$ and $\lVert\nabla_{y}g(\mathdutchcal{h}(u^{*}_{t},d_{t}))\rVert_{2}$ are justified, since $g(\cdot)$ is differentiable in a compact set. The persistency of excitation ensures that $\lVert\Delta u_{t}\rVert_{2}>0$ with high probability. Then, $\Sigma_{p,t}=\Sigma_{p_{1}}+\Sigma_{p_{2}}\lVert\Delta u_{t}\rVert_{2}^{2}\succ 0,\Sigma_{m,t}=\Sigma_{m_{1}}+\Sigma_{m_{2}}\lVert\Delta u_{t}\rVert_{2}^{2}+\Sigma_{m_{3}}\lVert\Delta u_{t}\rVert_{2}^{4}\succ 0$ if at least one $\Sigma_{p_{i}}\succ 0$ and one $\Sigma_{m_{j}}\succ 0$ for some $i,j$. Finally, even though the true sensitivity is state dependent, i.e., $H_{t}=\nabla_{u}\mathdutchcal{h}(u_{t},d_{t})$, the process and measurement noises in (6) and (8) allow to overapproximate the actual behavior of the sensitivity via these simplifications. In conclusion, Assumption 1 is reasonable. Then, with a persistently exciting $\Delta u$ as in (10), we have the following convergence result:

Proposition 1 establishes first that the estimated sensitivity $\hat{h}_{t}$ converges in expectation to the true sensitivity $h_{t}$ with a bounded covariance. Additionally, the control input $u_{t}$ converges to the AC-OPF solution $u_{t}^{*}$ from (2) with a quantifiable tracking error determined by the bound $C_{h,3}$ of the sensitivity estimation covariance, the variance $\sigma_{u}$ of the persistency of excitation noise $\omega_{u}$, and the temporal variation of the AC-OPF solution $\mathbb{E}[\lVert\Delta u^{*}_{t}\rVert_{2}^{2}]$, where $\Delta u^{*}_{t}$ can also be bounded by the temporal variation of $d_{t}$ and $\mathcal{U}_{t}$ in the AC-OPF (2) [28].

• •

  Distribution grid: We use the 3-phase, unbalanced IEEE 123-bus test feeder [19] in Figure 2.

• •

  Disturbance $d$: We consider uncontrollable active and reactive loads in our disturbance vector $d$. To generate these load profiles we use $1$-second resolution data of the ECO data set [29], then aggregate households and rescale them to the base loads of the 123-bus feeder. This gives us values of $d_{t}$ for every second during simulation time of 1h.

• •

  Controllable input set-points $u$: We add two solar PV systems and two wind turbines to the grid as in [8], see Figure 2. They can inject active power, and inject and absorb reactive power on all three phases, which gives us 24 control inputs. We consider a slack bus 150 in Figure 2, with a controllable voltage magnitude through, e.g., a tap changer, which makes in total $n_{u}=25$. The solar and wind generation profiles are generated based on a $1$-minute solar irradiation profile [30] and a $2$-minute wind speed profile [31]. Generation is assumed constant between samples. We use these profiles to set the time-varying upper limit of the feasible set $\overline{u}_{t}$, set the lower limit of active generation to $\underline{u}_{t}=0$, and define $\mathcal{U}_{t}=\{u\,|\,\underline{u}_{t}\leq u\leq\overline{u}_{t}\}$.

• •

  Output $y$: We consider as output $y$ the voltage magnitudes of all phases at all buses except the slack bus, given that it is a control input.

• •

  AC-OPF cost function in (2): We use a quadratic cost that penalizes deviating from a reference: $f(u)=\frac{1}{2}\lVert u-u_{\text{ref}}\rVert_{2}^{2}$. The reference $u_{\text{ref}}$ for the voltage magnitude at the slack bus is $1$ p.u. The reference for the controllable generation is the maximum installed power to promote using as much renewable energy as possible. The reference for reactive power is $0$. Note that the cost function is continuously differentiable, and has a strongly monotone and Lipschitz continuous gradient as required in Assumption 1. We consider the voltage limits $[0.94\,\text{p.u.},1.06\,\text{p.u.}]$ for all nodes as in [5, 8], and use the penalty function $g(y)=\frac{\rho}{2}\max\big{(}\left[\begin{smallmatrix}\mathbbm{1}\\ -\mathbbm{1}\end{smallmatrix}\right]y+\left[\begin{smallmatrix}-1.06\\ 0.94\end{smallmatrix}\right],0\big{)}^{2}$, with a sufficiently large penalization parameter $\rho=100$ to discourage violations. Again, this function is continuously differentiable, and has a monotone and Lipschitz continuous gradient.

• •

  Sensitivity process and measurement noises in (6) and (8): Under fast sampling rates $\Delta d_{t}$ may be negligible, especially when compared to $\Delta u_{t}$. Hence, for the simulation we assign $\Sigma_{p_{1}},\Sigma_{m_{1}},\Sigma_{m_{2}}$ to 0, and keep $\Sigma_{p_{2}},\Sigma_{m_{3}}\succ 0$. This ensures that $\Sigma_{p,t},\Sigma_{m,t}\succ 0$ for all $t$, as required by Assumption 1.

• •

  Persistency of excitation: We use a symmetric truncated Gaussian distribution with $\sigma_{u}=0.0001$ p.u. to introduce a low persistency of excitation noise $\omega_{u,t}$ that facilitates our sensitivity learning, but avoids introducing a big deviation in the input convergence, see (12).

• •

  Initializing sensitivity and linear model (3): We use the zero-injection operating point $u_{\text{op}}=0,d_{\text{op}}=0$ to initialize the sensitivity estimation, i.e., $\hat{H}_{0}=H_{0}=\nabla_{u}\mathdutchcal{h}(u,d)|_{(0,0)}$, see (3). In the first simulation (1: true admittance) we use the true admittance to compute $H_{0}$, in the second (2: perturbed admittance) we use a perturbed admittance matrix, where we have introduce an up to $20\%$ error in the admittance of the lines indicated in Figure 2.

We analyze the simulation performance of OFO with sensitivity learning (9) and (10), and compare it against an OFO with constant sensitivity (5). We validate both results in Proposition 1: First, the estimated sensitivity $\hat{H}_{t}$ converges to the real time-varying sensitivity $H_{t}$. Second, the input $u_{t}$ converges to the AC-OPF solution $u^{*}_{t}$ (2).