Tech note: smooth droop encoding

The Volt-var (VVC) and Volt-watt (VWC) characteristics of AS/NZS 4777.2 / IEEE 1547 are piecewise-linear (PWL) functions of voltage magnitude. PWL functions are not differentiable at their breakpoints, and primal–dual interior-point solvers such as Ipopt inherit Newton's requirement that the constraint functions be twice continuously differentiable: a kink makes the Jacobian jump and the Hessian a Dirac spike, which degrades the line search near a breakpoint. This note documents how BMOPFTools encodes the droop so that the OPF stays smooth, exact to a controllable tolerance, and numerically robust — and writes out the derivatives and the stable evaluation that the source only summarises.

The approach follows Mhanna, Geth, Quiertant & Mancarella [1], §II-C.1; the smoothing is the Chen–Harker–Kanzow–Smale (CHKS) softplus [2]. Numbered citations refer to the references below. The implementation lives in control_curves.jl.

1. The ReLU-sum encoding

Any continuous PWL function $f$ defined by non-decreasing breakpoints $(\bar x_1,\bar y_1),\dots,(\bar x_n,\bar y_n)$ and clamped flat outside $[\bar x_1,\bar x_n]$ can be written as a single baseline plus a sum of shifted/scaled rectified linear units $\operatorname{ReLU}(z)=\max(0,z)$:

\[f(U) \;=\; \bar y_1 \;+\; \sum_{i} a_i \,\operatorname{ReLU}(U-\bar x_i).\]

Each interior segment $i$ (between $\bar x_i$ and $\bar x_{i+1}$) with slope $s_i$ contributes two terms — $(+s_i,\bar x_i)$ turns the slope on at the segment start and $(-s_i,\bar x_{i+1})$ turns it off at the end — so the slopes telescope: the running slope below $\bar x_1$ is $0$, equals $s_i$ on segment $i$, and returns to $0$ above $\bar x_n$. This is the "two-triple per segment" form produced by breakpoints_to_triples.

For the canonical Volt-watt curve (100 % below 253 V dropping to 20 % at 260 V) the encoding is the two triples $(-\tfrac{0.8}{7},253)$ and $(+\tfrac{0.8}{7},260)$ on a baseline of $1.0$. Zero-slope (deadband) segments contribute no triples.

The encoding is exact for the kinked curve and is used as the reference (curve_value_exact) against which the smooth surrogate is tested.

2. The softplus surrogate

The non-smoothness is confined to the ReLU kink at $z=0$. Replace each ReLU by its softplus surrogate with smoothing parameter $\varepsilon>0$,

\[\operatorname{ReLU}^{\varepsilon}(z) \;=\; \varepsilon\,\log\!\bigl(1+e^{z/\varepsilon}\bigr),\]

giving the smooth droop

\[f^{\varepsilon}(U) \;=\; \bar y_1 \;+\; \sum_i a_i\,\varepsilon\,\log\!\bigl(1+e^{(U-\bar x_i)/\varepsilon}\bigr).\]

This function $f^{\varepsilon}$ is infinitely differentiable, $C^\infty$, and its slopes match the original droop exactly in the limit. Two properties make it a principled choice:

• Monotone, bounded error. The softplus brackets the ReLU from above with a uniform bound [1, 2]:

  \[0 \;\le\; \operatorname{ReLU}^{\varepsilon}(z) - \operatorname{ReLU}(z) \;\le\; \varepsilon\log 2 .\]

  The worst-case error is $\varepsilon\log 2$, attained at the kink, and it shrinks linearly in $\varepsilon$. Summing over triples, the curve error is bounded by $\bigl(\sum_i |a_i|\bigr)\varepsilon\log 2$.
• Exact limit. $\operatorname{ReLU}^{\varepsilon}(z)\to\operatorname{ReLU}(z)$ pointwise as $\varepsilon\to 0^+$, so the smoothing is a controllable approximation.

2.1. The native logistic / Swish alternative

For solver backends that expose logistic as a native nonlinear primitive, BMOPFTools also supports softplus=:swish. It replaces each hinge with

\[\operatorname{ReLU}^{\varepsilon}_{\mathrm{swish}}(z) = z\,\sigma\!\left(\frac{z}{\varepsilon}\right).\]

This is the same expression obtained by blending the two adjacent affine pieces with a logistic partition of unity. The breakpoint-to-triple encoding, relative ε scaling, and per-context operator cache are unchanged; only the smooth hinge operator changes. piecewise_linear_value selects the matching numeric oracle with encoding=:swish.

Swish is not a conservative replacement for softplus, and the two err in opposite directions. Softplus overestimates the hinge everywhere, with a one-sided error peaking at $\varepsilon\log 2 \approx 0.6931\varepsilon$ at $z = 0$. Swish underestimates it everywhere — $z\,\sigma(z/\varepsilon) \le \max(z, 0)$ for every $z$ — with a signed error that reaches $-0.2785\varepsilon$ on both sides of the hinge, at $z \approx \pm 1.2785\varepsilon$:

So Swish is the tighter pointwise approximation of the two; what it gives up is structure, not accuracy. Its derivative dips below zero and its second derivative changes sign, so it is neither monotone nor convex near a hinge.

A uniform underestimate of each hinge does not make the assembled curve a uniform underestimate: where a hinge slope $a_i$ is negative, the term $a_i \cdot (\text{negative error})$ is positive. A signed hinge sum can therefore overshoot a nominal clamp just below a breakpoint, or leak a small nonzero value into a deadband. The standard softplus remains the default.

The native mode emits a JuMP nonlinear expression whose operator head is :logistic. That is not one of MathOptInterface's default univariate operators, so the mode is strictly backend-specific: selecting softplus=:swish does not imply compatibility with Ipopt, DiffOpt, or any solver that does not advertise the logistic primitive — those reject the model with MOI.UnsupportedNonlinearOperator at optimize! rather than solving something else. Gurobi (12.0 or newer) is the motivating backend, where :logistic maps to GRB_OPCODE_LOGISTIC; it is also the only one of the three encodings Gurobi accepts, so see the solver guide before running a control-curve case there.

3. Closed-form first and second derivatives

Because the operator is a single-argument function, BMOPFTools registers it with JuMP's add_nonlinear_operator together with explicit first and second derivatives, avoiding automatic-differentiation overhead and giving Ipopt an exact Hessian. With $z=U-\bar x$ and the logistic (sigmoid) $\sigma(t)=\tfrac{1}{1+e^{-t}}$,

\[\frac{d}{dU}\operatorname{ReLU}^{\varepsilon}(z) = \sigma\!\Bigl(\tfrac{z}{\varepsilon}\Bigr), \qquad \frac{d^2}{dU^2}\operatorname{ReLU}^{\varepsilon}(z) = \frac{1}{\varepsilon}\,\sigma\!\Bigl(\tfrac{z}{\varepsilon}\Bigr)\Bigl(1-\sigma\!\Bigl(\tfrac{z}{\varepsilon}\Bigr)\Bigr).\]

The first derivative is a smooth switch rising from $0$ to $1$ across a band of width $\mathcal O(\varepsilon)$ around the breakpoint; the second derivative is a bump of height $\mathcal O(1/\varepsilon)$ and width $\mathcal O(\varepsilon)$ — a smooth, integrable stand-in for the Dirac curvature of the exact kink. For the full curve the derivatives are the corresponding $a_i$-weighted sums:

\[f^{\varepsilon\prime}(U) = \sum_i a_i\,\sigma\!\Bigl(\tfrac{U-\bar x_i}{\varepsilon}\Bigr), \qquad f^{\varepsilon\prime\prime}(U) = \frac{1}{\varepsilon}\sum_i a_i\,\sigma_i(1-\sigma_i),\quad \sigma_i=\sigma\!\Bigl(\tfrac{U-\bar x_i}{\varepsilon}\Bigr).\]

These are exactly the f, df, d2f closures passed to the operator in relu_operator.

4. Numerically stable evaluation (StatsFuns)

A literal evaluation of $\log(1+e^{t})$ and $\tfrac{1}{1+e^{-t}}$ with $t=z/\varepsilon$ overflows or underflows well inside the operating range — and the problem gets worse as $\varepsilon\to 0$, because $t=z/\varepsilon$ grows like $1/\varepsilon$:

• for large positive $t$, $e^{t}$ overflows to Inf, so $\log(1+e^t)$ returns Inf instead of $\approx t$;
• for large negative $t$, $e^{-t}$ overflows in the logistic denominator, and $1+e^{t}\to 1$ loses all the information in $e^t$.

BMOPFTools therefore evaluates the surrogate through StatsFuns.jl's log1pexp and logistic, which use the regime-split identities of Mächler [3]:

\[\log(1+e^{t}) = \begin{cases} e^{t} & t \lesssim -37 \quad(\text{1+}e^t\approx 1;\ \log 1p)\\ \log\!\bigl(1+e^{t}\bigr) & -37 \lesssim t \lesssim 18\\ t + e^{-t} & 18 \lesssim t \lesssim 33.3\\ t & t \gtrsim 33.3 \end{cases}\]

so the result is accurate to full Float64 precision and never overflows: the large-$t$ branch returns $t$ (the ReLU asymptote) directly, and the small-$t$ branch uses log1p to retain precision near zero. logistic applies the companion overflow-safe split, evaluating $\tfrac{e^{t}}{1+e^{t}}$ for $t<0$ and $\tfrac{1}{1+e^{-t}}$ for $t\ge 0$ so the exponential argument is always non-positive. The second derivative reuses $\sigma$, so it inherits the same stability. This is what lets the encoding use aggressively small $\varepsilon$ without the floating-point failures the naive formulas would suffer [1, §III-C].

5. Public API for custom control laws

The same encoding is available to staged-OPF extensions, so a downstream package can implement local control laws without copying BMOPFTools' nonlinear operator machinery. The numeric oracle evaluates either the exact controller or the same smooth approximation used in the optimisation model:

xs = [0.90, 0.98, 1.02, 1.10]
ys = [0.44, 0.0, 0.0, -0.60]

q_exact  = piecewise_linear_value(1.05, xs, ys)
q_smooth = piecewise_linear_value(1.05, xs, ys; epsilon=2e-3)
q_swish  = piecewise_linear_value(1.05, xs, ys;
                                  epsilon=2e-3, encoding=:swish)

Inside a model_hook!, pass a JuMP scalar expression and an absolute smoothing width in the model's working input units:

model_hook! = ctx -> begin
    model = opf_model(ctx)
    # `u` is a caller-constructed local voltage-magnitude expression in p.u.
    q_fraction = opf_piecewise_linear_expression(
        ctx, u, xs, ys; epsilon=2e-3)
    JuMP.@constraint(model, q == q_rating_pu * q_fraction)
end

All arguments to one curve are interpreted in working coordinates: input, breakpoints, and epsilon share input units, while values set the output units. Convert SI quantities with opf_bases when the staged model uses per-unit coordinates. Curve data are fixed finite numbers with strictly increasing breakpoints; changing their length or order requires rebuilding the JuMP graph.

Calls on the same context and epsilon share one cached smooth-ReLU operator. This keeps large populations of independently controlled devices compact and avoids re-registering identical nonlinear functions. The context's softplus mode is honoured automatically, including softplus=:builtin for DiffOpt workflows and softplus=:swish for native-logistic backends. The function only returns an expression: callers retain control over equality, inequality, current-limit, and topology-realizability constraints.

6. Choosing the smoothing $\varepsilon$

The smoothing is relative to the voltage scale. breakpoints_to_triples works in model units (per-unit when the OPF is solved per-unit), and the operator uses

\[\varepsilon = \texttt{relu\_eps}\times\overline{x},\qquad \overline{x}=\tfrac1n\sum_i \bar x_i,\]

i.e. relu_eps (default 2e-3, exposed as the volt_var_watt_eps keyword of solve_opf) times the mean breakpoint. Scaling $\varepsilon$ to the mean breakpoint keeps the corner radius a fixed fraction of the voltage range, so the droop is identical in SI and per-unit solves.

The tradeoff is the usual one:

• smaller $\varepsilon$ ⇒ tighter match to the standard's curve (error $\propto\varepsilon$), but a sharper second-derivative bump ($\propto 1/\varepsilon$) that can slow or trip the line search if a solution sits exactly on a breakpoint;
• larger $\varepsilon$ ⇒ rounder corners and easier conditioning, at the cost of a visibly smoothed deadband.

Mhanna et al. [1, §IV-C] benchmark this sweep and recommend a softplus smoothing small enough to stay well inside the $\pm 7.5\,\%$ voltage-accuracy band of real IBRs while preserving reliability; in their normalisation $\varepsilon\le 10^{-5}$ was the sweet spot. The default here is deliberately conservative and can be tightened per study. The same paper presents a quadratic Bézier spline as an alternative smoothing with an equivalent $\Delta U$ tolerance; BMOPFTools ships the softplus form, which has the cheaper gradient/Hessian evaluation.

Parameterized curve points

When a coefficient provider supplies breakpoints or ordinates, the model keeps the curve's nominal point count, hinge structure, and smoothing width fixed. This is essential for repeated parameter updates: the JuMP graph cannot grow or branch according to a parameter value. The engine retains every segment whose ordinate may change in the symbolic hinge sum. A segment with two fixed equal ordinates is omitted, avoiding a live zero-over-parameterized-gap expression. It also stamps strict ordering guards using a small gap derived from the nominal curve. Crossing two breakpoints therefore yields an infeasible solve rather than a different curve topology and an untrustworthy derivative.

For repeated coefficient updates, also follow the parameter-update workflow. Operator compatibility, evaluation stability, and cache freshness are separate concerns.

7. Provenance and related work

The encoding above is the composition of two well-established ideas; naming them clarifies what is standard and what is specific to this implementation.

Representing a PWL function as a sum of hinges. Writing a continuous PWL function as a baseline plus shifted $\max(0,\cdot)$ terms (§1) is the degree-1 truncated power basis of spline theory. It is the building block of Friedman's multivariate adaptive regression splines (MARS) [7] and of Breiman's hinging hyperplanes, and the canonical-PWL representation theory of Chua & Kang [8] shows that every continuous PWL function is a linear combination of maxima of affine functions — the single-variable case being exactly the hinge sum used here. Equivalently, a one-dimensional ReLU sum is a shallow ReLU network, which is why the encoding reads like one.

Smoothing the hinge with a softplus. Replacing $\max(0,z)$ by $\varepsilon\log(1+e^{z/\varepsilon})$ (§2) is the smoothing-of-the-plus-function idea of Chen & Mangasarian [9], who obtain the softplus by integrating a sigmoid precisely to convert nonsmooth inequality/complementarity problems into smooth ones solvable by Newton — the same motivation as here. The error bound and the broader smoothing family are the Chen–Harker–Kanzow–Smale (CHKS) line [2]. A unifying view: $\operatorname{softplus}(x)=\varepsilon\,\operatorname{logsumexp}(0,x/\varepsilon)$, and log-sum-exp is the canonical entropic (Nesterov) smoothing of $\max$ [10], so the surrogate is the simplest smooth-max. More recently, the pattern "relax to a program, then smooth it" has been systematised under differentiable programming [11], where softplus is a running example.

What is specific here. Neither ingredient is novel in isolation; the contribution is their exact-slope composition for standardised droop — the baseline and triples are placed analytically so the surrogate reproduces the AS/NZS 4777.2 / IEEE 1547 slopes in the $\varepsilon\to 0$ limit — together with the numerically stable log1pexp/logistic evaluation (§4) and deployment inside a neutral-explicit four-wire OPF [1]. The same smooth-droop strategy is applied to AC–DC converter droop with saturation in a companion paper [12].

References

1. S. Mhanna, F. Geth, L. Quiertant, P. Mancarella, "Volt-VAr-Watt Optimization in Four-Wire Low-Voltage Networks: Exact Nonlinear Models and Smooth Approximations," IEEE Trans. Power Systems, 2026, doi:10.1109/TPWRS.2026.3677246.
2. B. Chen, P. T. Harker, "Smooth approximations to nonlinear complementarity problems," SIAM J. Optim. 7 (2) (1997) 403–420.
3. M. Mächler, "Accurately computing $\log(1-e^{-|a|})$ — assessed by the Rmpfr package," 2012 (the log1pexp/logistic regime split used by StatsFuns.jl).
4. J. Huchette, J. P. Vielma, "Nonconvex piecewise linear functions: Advanced formulations and simple modeling tools," Operations Research 71 (5) (2023) 1835–1856 (the exact binary/logarithmic PWL alternative).
5. A. Wächter, L. T. Biegler, "On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming," Math. Program. 106 (1) (2006) 25–57 (Ipopt; the $C^2$ requirement).
6. M. Lubin, O. Dowson, J. Dias Garcia, J. Huchette, B. Legat, J. P. Vielma, "JuMP 1.0: Recent improvements to a modeling language for mathematical optimization," Math. Program. Comput., 2023 (user-defined nonlinear operators).
7. J. H. Friedman, "Multivariate adaptive regression splines," Ann. Statist. 19 (1) (1991) 1–67 (hinge-function / truncated-power basis).
8. L. O. Chua, S. M. Kang, "Section-wise piecewise-linear functions: Canonical representation, properties, and applications," Proc. IEEE 65 (6) (1977) 915–929 (canonical PWL representation).
9. C. Chen, O. L. Mangasarian, "A class of smoothing functions for nonlinear and mixed complementarity problems," Comput. Optim. Appl. 5 (2) (1996) 97–138 (softplus as the integral of a sigmoid, smoothing the plus function).
10. Yu. Nesterov, "Smooth minimization of non-smooth functions," Math. Program. 103 (1) (2005) 127–152 (entropic smoothing of $\max$ via log-sum-exp).
11. M. Blondel, V. Roulet, "The Elements of Differentiable Programming," arXiv:2403.14606, 2024 (smoothing programs; softplus as a running example).
12. G. Mohy-ud-din, R. Heidari, F. Geth, H. Ergun, S. M. M. Uddin, "AC-DC Power Systems Optimization with Droop Control Smooth Approximation," arXiv:2409.18376, 2024 (companion smooth-droop application).