[question] AC sensitivity analysis

[leuchtum]

@rdzman In this comment you shared the makePTDFac.m function. First of all, thanks for that, it really is what I'm looking for! But I have some questions about it:

1. Why isn't the function makePTDFac part of the official distribution? I came across this issue by pure luck as I was looking for ways to analyze sensitivity. Do you have any reservations about the function?
2. Im a bit confused about the resulting two (block)matrices. Is the following interpretation of the first one correct? (Ignoring slack just for simplicity)

H^f 
= B^f (B^{bus})^{\ -1}
= B^f J^{\ -1}
=
\left[ 
\begin{array}{c|c} 
  \frac{\delta P_{ij}}{\delta P_{i}} & \frac{\delta P_{ij}}{\delta Q_{i}} \\ 
  \hline 
  \frac{\delta Q_{ij}}{\delta P_{i}} & \frac{\delta Q_{ij}}{\delta Q_{i}}
\end{array} 
\right] 

3. Do you have any detailed explanation of the calculation? The math behind it is not entirely clear to me.

[rdzman]

1. I whipped this together a number of years ago, but never found the time to adequately test it. It appears I also started writing a test for it (I have a t_makePTDFac.m) that attempts to check the results using numerical perturbation of PF and OPF solutions, but the tests are not currently passing. I don't know if it's an issue with the tests or the function.
2. That looks correct to me.
3. I'm afraid I don't, but the rough idea is to define ...

B_f
= \left[
\begin{array}{c|c}
 \frac{\delta P_{ij}}{\delta \Theta_{i}} & \frac{\delta P_{ij}}{\delta V_{i}} \\ 
 \hline
 \frac{\delta Q_{ij}}{\delta \Theta_{i}} & \frac{\delta Q_{ij}}{\delta V_{i}}
\end{array}
\right]

(without the columns that correspond to the reference bus and buses with fixed voltage magnitude) and

B_\mathrm{bus}^{-1}
= \left[
\begin{array}{c|c}
 \frac{\delta P_{i}}{\delta \Theta_{i}} & \frac{\delta P_{i}}{\delta V_{i}} \\ 
 \hline
 \frac{\delta Q_{i}}{\delta \Theta_{i}} & \frac{\delta Q_{i}}{\delta V_{i}}
\end{array}\right]^{\ -1}
= \left[
\begin{array}{c|c}
 \frac{\delta \Theta_{i}}{\delta P_{i}} & \frac{\delta V_{i}}{\delta P_{i}} \\ 
 \hline
 \frac{\delta \Theta_{i}}{\delta Q_{i}} & \frac{\delta V_{i}}{\delta Q_{i}}
\end{array}\right]

(without the rows that correspond to the reference bus and buses with fixed voltage magnitude, and without columns corresponding to slack buses).

Then

H_f  = B_f B_\mathrm{bus}^{-1}
= \left[
\begin{array}{c|c}
 \frac{\delta P_{ij}}{\delta \Theta_{i}} & \frac{\delta P_{ij}}{\delta V_{i}} \\ 
 \hline
 \frac{\delta Q_{ij}}{\delta \Theta_{i}} & \frac{\delta Q_{ij}}{\delta V_{i}}
\end{array}
\right]
\left[
\begin{array}{c|c}
 \frac{\delta \Theta_{i}}{\delta P_{i}} & \frac{\delta V_{i}}{\delta P_{i}} \\ 
 \hline
 \frac{\delta \Theta_{i}}{\delta Q_{i}} & \frac{\delta V_{i}}{\delta Q_{i}}
\end{array}\right]
=
\left[ 
\begin{array}{c|c} 
 \frac{\delta P_{ij}}{\delta P_{i}} & \frac{\delta P_{ij}}{\delta Q_{i}} \\ 
 \hline 
 \frac{\delta Q_{ij}}{\delta P_{i}} & \frac{\delta Q_{ij}}{\delta Q_{i}}
\end{array} 
\right] 

Sorry, I don't have more, but hopefully this is at least a little helpful.

[leuchtum]

Thank you for the detailed answer, I really appreciate it. It's clear to me that the slack bus is removed, but why is it important to also remove the buses where the tension is fixed?

My understanding of your outline is to calculate the AC PTDF only w.r.t. PQ buses and ignoring slack and PV buses. But a branch is sensitive to a change in real power injection at a PV bus nevertheless, right?

[rdzman]

Yes, branch flows change with a changes in real power injections at PV buses. But what the sensitivities are depends on what assumptions you make about how voltages vary as you perturb the injections. If you assume all voltages are free to vary, you get one set of sensitivities. If you assume some voltages are fixed, you get another set of sensitivities.

At least, I believe that was my rationale for doing it the way I did, and for allowing the user to specify which voltages were fixed. I would love for someone to look at it closely and confirm (or not) that it is correct. If it is correct, one should be able to verify it via numerical perturbations.

Since I don't have time to dig into this in more detail now, I've attached the t_makePTDFac.m.txt that I started working on, in case it's helpful for anyone digging further. Be sure to remove the .txt extension.

[leuchtum]

@rdzman I went down a rabbit hole, but finally I got it. The error itself is related to testing and is suuuuper simple. In my opinion the makePTDFac is implemented correctly. All your assumptions above make total sense to me.

Here is a snippet from your test:

%% column 2 of PTDF reproducible by perturbation?
mpc2 = mpc0;
mpc2.gen(2, PG) = mpc2.gen(2, PG) + dP * baseMVA;
[baseMVA, bus, gen, branch, success, et] = runpf(mpc2, mpopt);
Ff2 = [ branch(:, PF); branch(:, QF) ] / baseMVA;
Ft2 = [ branch(:, PT); branch(:, QT) ] / baseMVA;
t_is((Ff2 - Ff0) / dP, Hf(:, 2), 5, 'Hf(:, 2) matches numerical PF perturbation');
t_is((Ft2 - Ft0) / dP, Ht(:, 2), 5, 'Ht(:, 2) matches numerical PF perturbation');

The index for the PTDF matrices needs to be selected from GEN_BUS:

%% column 2 of PTDF reproducible by perturbation?
mpc2 = mpc0;
mpc2.gen(2, PG) = mpc2.gen(2, PG) + dP * baseMVA;
[baseMVA, bus, gen, branch, success, et] = runpf(mpc2, mpopt);
Ff2 = [ branch(:, PF); branch(:, QF) ] / baseMVA;
Ft2 = [ branch(:, PT); branch(:, QT) ] / baseMVA;
connected_bus = mpc2.gen(2, GEN_BUS);
t_is((Ff2 - Ff0) / dP, Hf(:, connected_bus), 5, 'Hf(:, ?) matches numerical PF perturbation');
t_is((Ft2 - Ft0) / dP, Ht(:, connected_bus), 5, 'Ht(:, ?) matches numerical PF perturbation');

With this, the test passes.

How to proceeded now? Tidy up the function and its test for a PR?

[rdzman]

Certainly. That sounds good to me. Thanks for figuring this out!

[leuchtum]

I will create a PR towards the end of March, as I don't have time before then.