Closes #1320.
Closes #1470.

Improved version of #1470.
Related to #1409.

Changes proposed in this Pull Request

• Make the approximation of the loss parabola independent of s_nom_max. To be able to do this a bound for the approximation error is needed, which has the added benefit of giving more fine-grained control of the approximation quality. Here an error bound based on absolute and relative error tolerances is used. This fixes Changing s_nom_max affects line losses #1320.

• Use secants instead of tangents for the approximation. That way, the approximation error is zero at the "corners" of the stepwise linear approximation. Since the optimum of the linear program tends to be located at a corner, this will hopefully reduce the effective error of the loss approximation

• Introduce lines losses close to 0. At the moment, flows smaller than s_nom_max/(2*n_tangents) are lossless.

Example

Here is a plot of how the approximation with secants for an exemplary line looks like with atol = 1 #MW and rtol = 0.1. I chose the parameters such that the numbers of loss constraints is relatively small.

The approximation could be improved by setting the endpoint of the last secant to s_nom_max. However, that would re-introduce a dependency on s_nom_max and should ideally only be done when it can be guaranteed that s_nom_max does not change between consecutive runs of the model.

Questions

• Are we happy with the overestimation of losses that may occur at s_nom_max? Of course a similar approximation with tangents could be constructed. However, that would underestimate the error at s_nom_max.
• Should we use equality constraints instead of inequality constraints? The plot in Improving transmission loss approximation for DC power flow #1409 suggests that there are many snapshots when the losses are bigger than they should be, i.e. the model uses them to destroy energy.

Derivation of formulas for relative and absolute error

Find a piecewise linear approximation $l(x)$ of the loss parabola $L(x) := rx^2$ on $[0, X]$ with:

1. $x_0 = 0$
2. $l(x_i) = L(x_i) = r x_i^2$, where $x_i$ are the endpoints of the linear segments

We want the endpoints to be on the parabola, since the optimal solution tends to be on the corners of a simplex, and hence that is where the approximation should be best.

(2) implies that our approximation consists of secants. Denote the endpoints of a secant as $a,b$ then

$l_{ab}(x) = \frac{L(b) - L(a)}{b - a }(x-a) + L(a) = r\frac{b^2 - a^2}{b - a }(x-a) + ra^2 = r(b+a)(x-a) + ra^2 = r((a+b)x - ab)$

The absolute error $AE(x)$ between $l_{ab}(x)$ and $L(x)$ is

$AE(x) = r((a+b)x - ab) - rx^2$

This is positive for $a < x < b$, since there the secant is above the parabola. To compute it's maximum, remember that the derivative has to be 0.

$0 = \frac{d}{dx} AE(x) = r(a+b)x - 2rx \Rightarrow x_{max} = \frac{a+b}{2}$

Hence the error is maximal at the midpoint. Assume we specify the absolute error tolerance as $\varepsilon$. From this we can compute the interval length $b-a$

$\varepsilon = AE(\frac{a+b}{2}) = r((a+b)\frac{a+b}{2} - ab) - r(\frac{a+b}{2})^2 = r( \frac{(a+b)^2}{2} - \frac{(a+b)^2}{4} - ab) = r(\frac{(a+b)^2}{4} - \frac{4ab}{4}) = \frac{r(a-b)^2}{4}$

which implies $b(a, \varepsilon) = 2\sqrt{\frac{\varepsilon}{r}} + a$.

Hence, to construct an approximation with absolute error at most $\varepsilon$, we choose the endpoints of the secants at $x_i = \dots, -4\sqrt{\frac{\varepsilon}{r}}, -2\sqrt{\frac{\varepsilon}{r}}, 0, 2\sqrt{\frac{\varepsilon}{r}}, 4\sqrt{\frac{\varepsilon}{r}}, \dots$

For small error tolerances and large line flows we might need many secants, which would translate to many additional constraint and slow down solving. So, in this case we will use an approximation based on the relative error instead.

The relative error is

$RE(x) = \frac{r((a+b)x - ab) - rx^2}{rx^2} = \frac{a+b}{x} - \frac{ab}{x^2} -1$

It's derivative is $\frac{2ab - (a+b)x}{x^3}$, which is zero at $x_{max}=\frac{2ab}{a+b}$.
Assuming the relative error tolerance is $\delta$ we get:

$\delta = RE(\frac{2ab}{a+b})= \dots = \frac{(b-a)^2}{4ab}$

We want to build a step rule $b(a)$ to construct a new end point $b$ from a given starting point $a$. Hence we rewrite the equation above as a polynomial in $b$ with constants $a$ and $\delta$:

$\delta = \frac{(b-a)^2}{4ab} \Leftrightarrow b^2 - a(2+4\delta)b + a^2 = 0$

This polynomial has the roots $b = a(1 + 2\delta \pm 2\sqrt{\delta(1+\delta)})$. It can be shown that the root with - term is smaller than $a$. Since we want $b>a$ we take the other root and get

$b(a,\delta) = a(1 + 2\delta + 2\sqrt{\delta(1+\delta)})$

For small $a$ this gets very impractically small. Therefore we use the step based on absolute error for small $a$, and the step based on relative error for large $a$. Combined we have

$x_{i+1}(x_{i},\varepsilon,\delta) := b(a,\varepsilon,\delta) = \max(2\sqrt{\frac{\varepsilon}{r}}+ a, 2(\delta + \sqrt{\delta +\delta^2})a + a)$

We compute new endpoints until $x_I > X$. Then $l(x) = \sum_{i=0}^I l_{x_{i}x_{i+1}}(x)$.

Checklist

• Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in docs.
• Unit tests for new features were added (if applicable).
• A note for the release notes docs/release-notes.md of the upcoming release is included.
• I consent to the release of this PR's code under the MIT license.