Transmission network effects of wind forecast deviations

Jonas Kersulis

January 13, 2017

Motivation

Texas has >17,000 MW of installed wind capacity. (source)

November 2016: record peak wind generation of 15,033 MW

Texas: where the people are

(source)

Texas: where the wind is

(source)

Texas: moving wind power to people

(source)

Research question

Wind forecasts are often unreliable.

Most wind fluctuations are harmless, but certain wind patterns violate system constraints.

What are the smallest (most likely) problematic fluctuations (forecast deviations)?

Use optimization to rank dangerous patterns by likelihood.

New insight into system vulnerability allows system operators to be prepared.

Modeling

Engineers often make assumptions.

  • Underlying physics --> accurate model
  • Various assumptions --> intermediate simpler models
  • All the assumptions --> minimal model

Consider this pattern in power flow and transmission line rating.

Power flow

Accurate: AC power flow

Current: $I_{ik} = y_{ik}(E_i-E_k),$ where

$\begin{align} E_i &= V_i\angle\theta_i \\ E_k &= V_k\angle\theta_k \\ y_{ik} &= g_{ik} + jb_{ik} \end{align}$

Using $S=EI^*$, derive

$\begin{align} P_{ik} &= V_iV_k(G_{ik}\cos\theta_{ik} + B_{ik}\sin\theta_{ik}) \\ Q_{ik} &= V_iV_k(G_{ik}\sin\theta_{ik} - B_{ik}\cos\theta_{ik}), \end{align}$

where $G_{ik} + jB_{ik}$ is the admittance matrix element corresponding to line $ik$, and $\theta_{ik}=\theta_i-\theta_k$.

Minimal: DC power flow

Begin with AC:
$\begin{align} P_{ik} &= V_iV_k(G_{ik}\cos\theta_{ik} + B_{ik}\sin\theta_{ik}) \\ Q_{ik} &= V_iV_k(G_{ik}\sin\theta_{ik} - B_{ik}\cos\theta_{ik}), \end{align}$
Voltages are typically close to 1.0:
$\begin{align} P_{ik} &= G_{ik}\cos\theta_{ik} + B_{ik}\sin\theta_{ik} \\ Q_{ik} &= G_{ik}\sin\theta_{ik} - B_{ik}\cos\theta_{ik}, \end{align}$
Reactive flow depends mainly on voltage difference, let's ignore it:
$\begin{align} P_{ik} &= G_{ik}\cos\theta_{ik} + B_{ik}\sin\theta_{ik} \end{align}$
Well, resistance is also pretty small...
$\begin{align} P_{ik} &= B_{ik}\sin\theta_{ik} \end{align}$
Know what else is small? Angle differences.
$\begin{align} P_{ik} &= B_{ik}\theta_{ik} \end{align}$

...and we've run out of assumptions to make.

Intermediate: AC pf relaxations

Lots of literature and active research.

Skipping this today.

Why DC power flow?

  • DC PF: Linearity means fast, guaranteed solutions!
  • AC PF: tougher to compute, may fail to converge at any time.

Power flow is often embedded in our algorithms. Unit commitment, economic dispatch, the algorithm I'll present later, etc.

  • DC PF: linear constraints are easy to deal with.
  • AC PF: nonlinear, non-convex constraints are not.

Line rating

Accurate: IEEE Standard 738

Heat balance equation: \begin{align} \frac{dT_\text{avg}}{dt} &= \frac{1}{mC_p}\left(I^2R(T_\text{avg}) - q_c - q_r + q_s\right) \end{align}

  • $T_{avg}$ is conductor average line temperature
  • $mC_p$ is heat capacity
  • $R(T_{avg})$ emphasizes dependence of resistance on temperature
  • convection heat loss $q_c$ varies with
    • conductor diameter
    • wind velocity
    • air density and viscosity
  • radiated heat loss $q_r$ varies with conductor emissivity and diameter according to Stefan-Boltzmann equation
  • solar heat gain $q_s$ varies with geometry, solar intensity, and line reflectivity

Minimal: active power rating

Assume maximum sun, minimum wind. Determine how much active power can flow indefinitely at some conservative power factor without overheating the line.

Dead simple, but leaves lots of headroom.

Intermediate 1: approximate line current

Begin with the minimal line rating model: constrain DC-approximate active power $P_{ik} = B_{ik}\theta_{ik}$.

From the heat balance equation, line temperature depends on $I_{ik}^2$, not just $P_{ik}$. $I_{ik}$ varies with voltage magnitudes $V_i$ and $V_k$ in addition to $\theta_{ik}$.

Line current approximation that is sensitive to voltage magnitudes:

\begin{align} \lvert I_{ik} \rvert^2 &\approxeq B_{ik} \left( V_i^2 + V_k^2 - 2V_iV_k\cos\theta_{ik} \right),\quad \theta_{ik} = \theta_i - \theta_k \end{align}

Intermediate 2: linearized heat balance equation

Begin with heat balance equation.

$\begin{align} \frac{dT_{\text{avg},ik}}{dt} &= \frac{1}{mC_p}\left(I^2R(T_\text{avg,ik}) - q_c - q_r + q_s\right) \end{align}$

Approximate $I^2R(T_{avg})$ (see Almassalkhi and Hiskens).

$\begin{align} \frac{dT_{\text{avg},ik}}{dt} &= \frac{1}{mC_p}\left(r_{ik}\left(\frac{\theta_{ik}}{x_{ik}}\right)^2 - q_c - q_r + q_s\right) \end{align}$

Only $q_r$ is nonlinear in temperature. Linearize about the temperature limit, and the heat balance equation is trivially integrable: $$T_{\text{avg},ik}(t) \approx ce^{at} + k$$

Model combinations

DC powerflow AC powerflow
Active power constraint Done by Chertkov et al., covered today Done, Baghsorkhi and Hiskens
Approximate current constraint Done, covered today Simple modification of above
Linearized heat balance (temporal) Done, covered today No known solution method, likely intractable
IEEE 738 standard Future work? Holy grail

Active power constraint with DC power flow

Consider a transmission network with a significant number of wind farms.

Suppose the demand forecast is accurate, and generation dispatch does not change.

Suppose the wind forecast is not entirely accurate. Let each wind farm's active power output be an optimization variable.

Optimization problem

  • Objective
  • System constraints
  • Strain constraint

Objective: minimize the normed distance between the vector of wind optimization variables and their forecast values.

\begin{align} \min & \frac{1}{2} \left( R - R^0 \right)^\top \Lambda \left( R - R^0 \right) \end{align}

This ensures we find a likely dangerous pattern. (Note the covarianace matrix $\Lambda$.)

System constraints: power balance, angle reference, generator response parameter $\alpha$.

\begin{align} \text{subject to:}& \\ \sum_k( Y_{ik} \theta_k) & = G_{i} + R_{i} - D_{i} && \forall i \in N \\ G_i &= G_{i}^0 + k_i\alpha && \forall i \in N_g \\ \theta_{ref} & = 0 && \text{ (angle reference)} \\ R_{i} &\geq 0 && \forall i \in N_r \\ R_{i} &= 0 && \forall i \in N\backslash N_r \end{align}

Strain constraint: We add one additional constraint to force a line to saturate.

\begin{align} \theta_i - \theta_k & = x_{ik} P_{ik}^{lim} && \text{for each } (i,k) \in G \\ \end{align}

Outcome

We obtain an instanton candidate (or "dangerous wind pattern") for each line in the network.

Each candidate has a "score" (its distance from the wind forecast).

We call the candidate with the lowest score the instanton.

What does this look like?

Suppose we have a network with three nodes and three edges.

Assume:

  • DC power flow approximation holds
  • The conventional generator takes up all slack and is the angle reference

Variables

  • Two wind output decision variables: $R_{1}$ and $R_{2}$
  • Two bus phase angle variables: $\theta_1$ and $\theta_2$

Power balance equations:

$$\begin{align} R_{1} - d_1 &= y_{12}(\theta_1 - \theta_2) + y_{1S}\theta_1 \\ R_{2} - d_2 &= y_{12}(\theta_2 - \theta_1) + y_{2S}\theta_2 \\ \end{align}$$

Line constraints:

$$\begin{align} |y_{12}(\theta_1 - \theta_2)| &\leq P^{lim} \implies & |\theta_1 - \theta_2| &\leq \frac{P^{lim}}{y_{12}} \\ \\ |y_{1S}(\theta_1)| &\leq P^{lim} \implies & |\theta_1| &\leq \frac{P^{lim}}{y_{1S}} \\ \\ |y_{2S}(\theta_2)| &\leq P^{lim} \implies & |\theta_2| &\leq \frac{P^{lim}}{y_{2S}} \end{align}$$

Feasible region of optimization

Let's plug in some numbers and plot this set of inequalities.

Plim = 1;
    y12 = 0.6;
    y1S = 0.3;
    y2S = 0.2;

    R1 = 0.3;
    R2 = 1.4;

Approximate current line rating with DC powerflow

Recall that the DC power flow approximation involves four assumptions:

  • No line resistance
  • Small angle differences
  • Flat voltage profile
  • No reactive power flows

We can improve solution accuracy by getting rid of the last two assumptions.

Approximate line current

DC power flow assumes a "flat" voltage profile, but we have voltage magnitude data.

Let's put it to use. From before:

\begin{align} \lvert I_{ik} \rvert^2 &\approxeq B_{ik} \left( V_i^2 + V_k^2 - 2V_iV_k\cos\theta_{ik} \right),\quad \theta_{ik} = \theta_i - \theta_k \end{align}

We still assume no resistance, but now we allow for non-flat voltage profiles.

How good is the approximation?

Let's plot the exact and approximate current expressions on the same chart.

$$\begin{align} I_{ik} &= Y_{ik}(E_i-E_k) &\qquad \lvert I_{ik} \rvert^2 &\approxeq B_{ik} \left( V_i^2 + V_k^2 - 2V_iV_k\cos\theta_{ik} \right) \end{align}$$

If resistance is zero:

Perfect overlap.

Now let's introduce non-zero resistance. A typical line has $\lvert B_{ik}/G_{ik} \rvert \approx 4$, so let's plot that case:

Our approximation is good as long as line resistance isn't too big.

Why use this approximation?

It would make no sense to use the exact expression. The rest of our analysis assumes $r=0$.

More importantly, our current flow approximation is better than the DC approximation.

Any line with a significant voltage magnitude difference across it is more vulnerable than linear (DC) instanton analysis tells us.

Implementation and results

Solve the approximate current expression for angle difference:

\begin{align} \lvert I_{ik}^2 \rvert &\approxeq \frac{V_i^2 + V_k^2 - 2V_iV_k\cos \theta_{ik}}{x_{ik}^2} \\ \implies \theta_i - \theta_k &= \cos^{-1}\left( \frac{V_i^2 + V_k^2 - (\lvert I_{ik}^{lim} \rvert x_{ik})^2}{2V_iV_k} \right) \end{align}

Same form as the DC instanton strain constraint:

\begin{align} \theta_i - \theta_k & = x_{ik} P_{ik}^{lim} \end{align}

Simply swap out this DC constraint for our approximate current one.

Compare: DC results

Compare: approx. $I$ results

Differences are due to the non-flat voltage profile of the RTS-96:

Linearized heat balance equation with DC power flow

3-word review of work we just covered, from John Adams of ERCOT:

"I don't care."

System operators allow temporary line overloads all the time.

They understand that lines trip due to heat-induced sag, not instantaneous current (IEEE Standard 738).

In an operating environment, we must

  • keep track of line temperature
  • Accommodate changes in power flow and ambient conditions

This motivates a new framework.

Instantaneous to temporal

Previous work: instantaneous

  • Minimize forecast deviations
  • Power balance
  • Single angle difference

New framework: temporal

  • Minimize forecast deviations across multiple time steps
  • Power balance
  • Line temperature constraint (heat added across multiple time steps)

What makes it harder?

Instantaneous analysis: Quadratic program with analytic solution.

Temporal analysis: Non-convex QCQP with multiple solutions.

Remainder of talk

  • Thermal constraint
  • Temporal instanton model
  • Temporal optimization problem
  • Solution
  • Results

Thermal constraint

Losses

Heat is a result of line losses.

Starting with the AC line loss expression, Almassalkhi and Hiskens (2014) derived: $$\begin{align} f_{ij}^{\text{loss}} &\approx r_{ij}\left(\frac{\theta_{ij}}{x_{ij}}\right)^2 \end{align}$$

Here we make three assumptions:

  1. Voltage magnitudes are all 1 pu
  2. Cosine is represented by its second-order Taylor expansion
  3. $x_{ij} \geq 4r_{ij}$ (line reactance is at least four times as great as resistance).

Skipping a long derivation...

We can constrain the line's temperature at time $t_{n_t}$ to be equal to this limiting value by enforcing: \begin{align} \sum_{k=1}^{n_t} \left(\hat{\theta}_{ij}(t_k)\right)^2 &= \frac{c_1}{c_4}\left(T_\text{lim} - c_7\right), \end{align} where \begin{align} \hat{\theta}_{ij}(t_{k}) &= \theta_{ij}(t_k)\sqrt{ (e^{c_1\bar{t}})^{n_t-k+1} - (e^{c_1\bar{t}})^{n_t-k} },\end{align} and $c_i$ are constants.

Optimization problem

  • Objective
  • System constraints
  • Strain constraint

Objective: minimize the normed distance between the vector of wind optimization variables and their forecast values.

\begin{align} \min & \sum_{t=1}^{n_t} dev_t^\top Q_{dev} dev_t \end{align}

(Note that $Q_{obj}$ may be identity or could encode correlation.)

System constraints: power balance, angle reference, generator response parameter $\alpha$.

\begin{align} \text{subject to:} \\ \sum_k Y_{ij} \theta_{ij,t} & = (G_{i,t} + \rho_{i,t} + x_{i,t}) - D_{i,t} ~ \forall i \in \mathcal{N},~t\in 1... n_t \\ G_t &= G_{0,t} + k\alpha_t \quad \forall t\in 1\ldots n_t \\ \theta_{ref,t} & = 0 \quad \forall t\in 1\ldots n_t \\ \end{align}

Strain constraint: Force a line to sag.

\begin{align} T_{ij}[n_t] &= T_{ij}^\text{lim}\quad \text{for some }(i,j)\in \mathcal{G} \end{align}

Strain constraint: Force a line to sag.

\begin{align} T_{ij}[n_t] &= T_{ij}^\text{lim}\quad \text{for some }(i,j)\in \mathcal{G} \end{align}

From what we just discussed, this is really a norm constraint on a vector of scaled angle differences:

\begin{align*} \lVert \hat{\theta}_{ij}\rVert^2 &= c^2 \end{align*}

Recast as optimization problem

General form: quadratic objective function with linear and quadratic constraints. Stacking all variables into a single vector $z$, we find:

\begin{align} && \min z^\top Q_{obj} z &\\ & s.t. & Az &= b \\ && z^\top Q_{c}z &= c \end{align}

Let's consider each component:

  • Variables: $z$
  • Objective: $Q_{obj}$
  • Linear constraints: $A$ and $b$
  • Quadratic constraint: $Q_c$

Variables: $z$

Stack variables from all time steps. Put $\hat{\theta}_{ik}$ at the end. \begin{align} z &= \begin{bmatrix} dev^{(1)\top} & \theta^{(1)\top} & \alpha^{(1)} & \cdots & dev^{(n_t)\top} & \theta^{(n_t)\top} & \alpha^{(n_t)} & \hat{\theta}_{ik}^\top \end{bmatrix}^\top \end{align}

Objective: $Q_{obj}$

The objective weights only deviation variables $dev$.

Example: suppose there are two time steps and correlation is represented by $Q_{dev}$. Then

$$ z = \begin{bmatrix} dev^{(1)\top} & \theta^{(1)\top} & \alpha^{(1)} & dev^{(2)\top} & \theta^{(2)\top} & \alpha^{(2)} & \hat{\theta}_{ik}^\top \end{bmatrix}^\top $$

and $Q_{obj}$ is given by:

\begin{align*} Q_{obj} = \begin{bmatrix} Q_{dev} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & Q_{dev} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \end{align*}

Linear constraints: $A$ and $b$

$A$ has a block diagonal structure where each block is $(n_b+1)\times(n_{d}+n_b)$.

  • The first $n_b$ rows describe power balance and distributed slack behavior. We have the following power balance equation for each node:
\begin{align*} -dev_{i,t} + \sum\limits_{j} Y_{ij}\theta_{j,t} - k_i\alpha_{i,t} &= G_{i,t}^0 + \rho_{i,t} - D_{i,t} \end{align*}
  • The last row in each group of $(n_b+1)$ fixes the angle reference to zero.

Quadratic constraint: $Q_c$

Recall that the line temperature constraint in terms of $\hat{\theta}_{ij}$ is \begin{align*} \lVert \hat{\theta}_{ij}\rVert^2 &= c^2~. \end{align*}

Because the last $n_t$ elements of $z$ contain $\hat{\theta}_{ij}$, $Q_c$ is all zeros with a $n_t\times n_t$ identity matrix in the lower-right corner.

\begin{align} Q_c &= \begin{bmatrix} 0 & 0 \\ 0 & I\end{bmatrix} \end{align}

Solution

  • Quadratic equality constraint makes the problem non-convex.
  • It does look nice though.

Our QCQP is similar to the "trust-region subproblem". Dr. Dan Bienstock's insight from this related problem helped us derive a solution for ours.

Here's how we solve it:

  1. Translation to change $Az=b$ into $Ay=0$
  2. Kernel mapping to render linear constraints implicit
  3. Rotation to restore constraint's norm structure
  4. Elimination of unconstrained variables
  5. Solution via enumeration

Skipping a lot here! Paper coming soon...

Results for wind-augmented RTS-96 network

Acknowledgements

  • Dr. Ian Hiskens
  • Dr. Misha Chertkov
  • Dr. Scott Backhaus
  • Dr. Dan Bienstock

References

[1] M. Chertkov, F. Pan, and M. Stepanov, “Predicting failures in power grids: The case of static overloads,” IEEE Transactions on Smart Grid, vol. 2, no. 1, pp. 162–172, Mar. 2011.

[2] M. Chertkov, M. Stepanov, F. Pan, and R. Baldick, “Exact and efficient algorithm to discover extreme stochastic events in wind generation over transmission power grids,” in Proc. 2011 50th IEEE Conference on Decision and Control and European Control Conference (CDC-ECC), 2011, pp. 2174–2180.

[3] S. Baghsorkhi and I. Hiskens, “Analysis tools for assessing the impact of wind power on weak grids,” in Proc. Systems Conference (SysCon), 2012 IEEE International, 2012, pp. 1–8.

[4] H. Banakar, N. Alguacil, and F. Galiana, “Electrothermal coordination part I: theory and implementation schemes,” IEEE Transactions on Power Systems, vol. 20, no. 2, pp. 798–805, May 2005.

[5] “IEEE standard for calculating the current-temperature of bare overhead conductors,” IEEE Std 738-2012

[6] M. Almassalkhi and I. Hiskens, “Model-predictive cascade mitigation in electric power systems with storage and renewables – part I: Theory and implementation,” IEEE Transactions on Power Systems, vol. PP, no. 99, pp. 1–11, 2014.

[7] O. Mehanna, K. Huang, B. Gopalakrishnan, A. Konar, and N. Sidiropou- los, “Feasible point pursuit and successive approximation of non-convex QCQPs,” IEEE Signal Processing Letters, vol. PP, no. 99, pp. 1–1, 2014.

[8] D. Bienstock and A. Michalka, “Polynomial Solvability of Variants of the Trust-region Subproblem,” in Proceedings of the Twenty-Fifth Annual ACM-SIAM Symposium on Discrete Algorithms, ser. SODA ’14. Portland, Oregon: SIAM, 2014, pp. 380–390. [Online]. Available: http://dl.acm.org/citation.cfm?id=2634074.2634102

[9] H. Pandzic, Y. Dvorkin, T. Qiu, Y. Wang, and D. Kirschen. Unit Commitment under Uncertainty - GAMS Models, Library of the Renewable Energy Analysis Lab (REAL), University of Washington, Seattle, USA. [Online]. Available: http://www.ee.washington.edu/research/real/gams_code.html