Grey-Box Design: Sizing a Pipeline Around an Iterative Friction Solver#
External Functions covers how dm.external works:
its shape rule, how the derivatives get through, and its refusals. This notebook
uses it on an engineering problem, with a grey box of the kind that actually turns
up in practice. The block’s value comes out of a numerical solve, and its
derivatives come from implicit differentiation of that solve, not from any
closed-form expression.
The problem: choose the diameters of a three-segment water pipeline to minimize
annualized cost, trading pipe capital (which grows with diameter) against pumping
power (which falls steeply with it), subject to a pressure-drop budget along the
line. The friction factor in every segment is given by the Colebrook equation
[Colebrook, 1939], which is implicit in the friction factor and is
normally solved iteratively. That solver is the grey box; the rest of the model
is ordinary algebra [Biegler, 2010]. This is the same split Pyomo’s
ExternalGreyBoxModel handles [Hart et al., 2017].
The grey box: an iterative Colebrook solver#
For a pipe of diameter \(D\) and roughness \(\varepsilon\) carrying a flow at Reynolds number \(\mathrm{Re}\), the Darcy friction factor \(f\) satisfies
Write \(s = 1/\sqrt f\). With the volumetric flow \(Q\) fixed per segment, \(\mathrm{Re} = 4\rho Q/(\pi\mu D)\), so \(2.51/\mathrm{Re} = cD\) with \(c = 2.51\pi\mu/(4\rho Q)\), and the equation becomes the residual
The function below solves \(G = 0\) with Newton’s method. It uses a Python loop with a data-dependent exit, so it cannot be traced, and in a real code it might just as well be a call into a compiled property library.
import itertools
import time
import warnings
import discopt.modeling as dm
import numpy as np
warnings.filterwarnings("ignore")
# Water at 20 C in commercial steel pipe.
RHO, MU, ROUGH = 1000.0, 1.0e-3, 4.5e-5 # kg/m^3, Pa s, m
L = np.array([400.0, 250.0, 300.0]) # segment lengths, m
Q = np.array([0.05, 0.03, 0.02]) # segment flows, m^3/s (offtakes between segments)
K = 2.0 / np.log(10.0)
E = ROUGH / 3.7
C = 2.51 * np.pi * MU / (4.0 * RHO * Q) # 2.51/Re = C*D
def solve_colebrook(D, c, tol=1e-14, max_iter=50):
# Newton on G(s) = s + k ln(e/D + c D s); returns s = 1/sqrt(f).
s = 7.0
for _ in range(max_iter):
u = E / D + c * D * s
step = (s + K * np.log(u)) / (1.0 + K * c * D / u)
s -= step
if abs(step) < tol * max(1.0, abs(s)):
return s
raise RuntimeError(f"Colebrook did not converge at D={D}")
D_probe = np.array([0.20, 0.15, 0.12])
Re = 4 * RHO * Q / (np.pi * MU * D_probe)
f = np.array([solve_colebrook(d, c) ** -2 for d, c in zip(D_probe, C)])
print("Re =", np.round(Re))
print("f =", np.round(f, 5))
Re = [318310. 254648. 212207.]
f = [0.01634 0.01725 0.01804]
Derivatives by the implicit function theorem#
To use the solver inside an NLP we need \(df/dD\) and \(d^2f/dD^2\). Differentiating through the Newton iterations would be wrong (it differentiates the iteration count and the stopping test) and wasteful. The implicit function theorem differentiates the converged residual instead [Griewank and Walther, 2008]. From \(G(s(D), D) = 0\):
and then, since \(f = s^{-2}\),
The partials of \(G\) are elementary:
\(G_s = 1 + kcD/u\)
\(G_D = k\,u_D/u\), where \(u_D = -e/D^2 + cs\)
\(G_{ss} = -k(cD)^2/u^2\)
\(G_{Ds} = kc\,(u - Du_D)/u^2\)
\(G_{DD} = k\,(2eu/D^3 - u_D^2)/u^2\)
Each segment’s friction factor depends only on its own diameter. The Jacobian of
the block’s three outputs is therefore diagonal, and the Hessian has only its
\([i,i,i]\) entries nonzero. The shapes still follow dm.external’s rule:
jac is (3,) + (3,) and hess is (3,) + (3,) + (3,).
The counters let us see later how often the solver calls the block.
CALLS = {"fn": 0, "jac": 0, "hess": 0}
def _friction_and_derivatives(D, c):
s = solve_colebrook(D, c)
u = E / D + c * D * s
uD = -E / D**2 + c * s
Gs = 1.0 + K * c * D / u
GD = K * uD / u
Gss = -K * (c * D) ** 2 / u**2
GDs = K * c * (u - D * uD) / u**2
GDD = K * (2.0 * E * u / D**3 - uD**2) / u**2
s1 = -GD / Gs
s2 = -(GDD + 2.0 * GDs * s1 + Gss * s1**2) / Gs
return s**-2, -2.0 * s**-3 * s1, 6.0 * s**-4 * s1**2 - 2.0 * s**-3 * s2
def friction(D):
CALLS["fn"] += 1
return np.array([_friction_and_derivatives(float(d), c)[0] for d, c in zip(D, C)])
def friction_jac(D):
CALLS["jac"] += 1
return np.diag([_friction_and_derivatives(float(d), c)[1] for d, c in zip(D, C)])
def friction_hess(D):
CALLS["hess"] += 1
H = np.zeros((3, 3, 3))
for i, (d, c) in enumerate(zip(D, C)):
H[i, i, i] = _friction_and_derivatives(float(d), c)[2]
return H
Check the derivatives before trusting them#
discopt checks the shapes your callables return, but it has no way to check their values: a wrong Jacobian with the right shape just sends the NLP solver the wrong way. Compare against central differences before building any model. This is cheap, and it is the single most useful habit with grey boxes.
h = 1e-6
I3 = np.eye(3)
J_fd = np.column_stack(
[(friction(D_probe + h * e) - friction(D_probe - h * e)) / (2 * h) for e in I3]
)
H_fd = np.stack(
[(friction_jac(D_probe + h * e) - friction_jac(D_probe - h * e)) / (2 * h) for e in I3], axis=-1
)
J, H = friction_jac(D_probe), friction_hess(D_probe)
err_J = np.abs(J_fd - J).max() / np.abs(J).max()
err_H = np.abs(H_fd - H).max() / np.abs(H).max()
print(f"jac shape {J.shape} max rel. error vs FD: {err_J:.1e}")
print(f"hess shape {H.shape} max rel. error vs FD: {err_H:.1e}")
assert np.allclose(J_fd, J, rtol=1e-6) and np.allclose(H_fd, H, rtol=1e-5, atol=1e-8)
jac shape (3, 3) max rel. error vs FD: 1.8e-09
hess shape (3, 3, 3) max rel. error vs FD: 1.1e-10
The design model#
The pressure drop in segment \(i\) is given by Darcy–Weisbach with \(v = 4Q/(\pi D^2)\):
The annualized cost is the pipe capital, \(a\,L_iD_i^{1.5}\), plus electricity for the pump, \(b\,Q_i\Delta p_i\). Here \(b\) prices one watt of hydraulic power over 8000 h/yr at $0.12/kWh with a 75%-efficient pump. The whole line has a 3 bar budget, \(\sum_i \Delta p_i \le 3\times10^5\) Pa, which couples the three segments.
Only \(f\) is opaque. The \(D^{-5}\) law, the costs and the budget are written in dm.*,
so the rest of the model stays visible to discopt.
A_CAP = 1200.0 # $/yr per m of pipe per m^1.5 of diameter
B_POW = 8000 * 0.12 / 0.75 / 1000 # $/yr per W of hydraulic power
DP_MAX = 3.0e5 # Pa
G_DP = 8 * RHO * L * Q**2 / np.pi**2 # dp_i = f_i * G_DP[i] / D_i^5
colebrook = dm.external(
friction, jac=friction_jac, hess=friction_hess, shape=(3,), name="colebrook"
)
def build_pipeline(friction_factors):
# friction_factors(m, D) returns the three friction-factor expressions.
m = dm.Model("pipeline")
D = m.continuous("D", shape=(3,), lb=0.08, ub=0.5)
f = friction_factors(m, D)
dp = [f[i] * G_DP[i] / D[i] ** 5 for i in range(3)]
m.minimize(sum(A_CAP * L[i] * D[i] ** 1.5 + B_POW * Q[i] * dp[i] for i in range(3)))
m.subject_to(sum(dp) <= DP_MAX, name="pressure_budget")
return m, D
for k in CALLS:
CALLS[k] = 0
m_ext, D_ext = build_pipeline(lambda m, D: colebrook(D))
t0 = time.perf_counter()
res_ext = m_ext.solve(time_limit=120)
t_ext = time.perf_counter() - t0
D_opt = np.asarray(res_ext.value(D_ext), dtype=float)
dp_opt = friction(D_opt) * G_DP / D_opt**5
print(f"status : {res_ext.status}")
print(f"cost : ${res_ext.objective:,.2f}/yr ({t_ext:.1f} s)")
print(f"D : {np.round(D_opt, 4)} m")
print(f"dp : {np.round(dp_opt / 1e5, 3)} bar, total {dp_opt.sum() / 1e5:.3f} bar")
print(f"calls : {CALLS}")
Rust model repr unavailable — FBBT/presolve disabled: TypeError: Unknown expression type: CustomCall
No valid dual bound was produced for model 'pipeline': the relaxation layer could not bound this objective (status=feasible). The result cannot be certified globally optimal. A nonlinear term with no envelope is the usual cause; an epigraph reformulation (minimize z subject to f(x) <= z) often gives the relaxation something to bound.
status : feasible
cost : $72,349.67/yr (0.8 s)
D : [0.1621 0.1314 0.1118] m
dp : [1.189 0.806 1.005] bar, total 3.000 bar
calls : {'fn': 13, 'jac': 12, 'hess': 12}
The budget is active (total drop 3.000 bar). The counters show that discopt evaluated the friction solver, its Jacobian and its Hessian, which confirms that your derivative code is what the NLP used.
As with every dm.external model, the result is feasible, not optimal. The
NLP solver [Wächter and Biegler, 2006] found a point satisfying the KKT conditions, but
discopt has no relaxation of an opaque block, so it reports no bound:
print(f"bound={res_ext.bound} gap={res_ext.gap} gap_certified={res_ext.gap_certified}")
bound=None gap=None gap_certified=False
Was that the global optimum? Ask a symbolic twin#
For a real simulator you could not answer that. Colebrook is unusual in that it can be written algebraically, by making \(s_i\) a variable and imposing \(G(s_i, D_i) = 0\) as a constraint. discopt can then relax and branch on the whole model and certify a global optimum. That gives us an independent check on the grey-box answer.
def symbolic_friction(m, D):
s = m.continuous("s", shape=(3,), lb=2.0, ub=20.0)
for i in range(3):
m.subject_to(s[i] + K * dm.log(E / D[i] + C[i] * D[i] * s[i]) == 0)
return [1.0 / s[i] ** 2 for i in range(3)]
m_sym, D_sym = build_pipeline(symbolic_friction)
t0 = time.perf_counter()
res_sym = m_sym.solve(time_limit=300)
t_sym = time.perf_counter() - t0
D_twin = np.asarray(res_sym.value(D_sym), dtype=float)
print(
f"symbolic twin : {res_sym.status:8s} cost ${res_sym.objective:,.4f} "
f"gap_certified={res_sym.gap_certified} ({t_sym:.1f} s)"
)
print(
f"grey box : {res_ext.status:8s} cost ${res_ext.objective:,.4f} "
f"gap_certified={res_ext.gap_certified} ({t_ext:.1f} s)"
)
rel_diff = abs(res_ext.objective - res_sym.objective) / res_sym.objective
print(f"relative objective difference : {rel_diff:.1e}")
print(f"max diameter difference : {np.abs(D_opt - D_twin).max():.1e} m")
symbolic twin : optimal cost $72,349.6670 gap_certified=True (49.9 s)
grey box : feasible cost $72,349.6670 gap_certified=False (0.8 s)
relative objective difference : 1.2e-11
max diameter difference : 6.3e-12 m
The two agree to solver tolerance, and the twin’s answer carries a global certificate. The grey-box answer was globally optimal, but we know that only because of the twin. The certificate belongs to the twin and covers this instance only. On a different flow set or cost model the local solve could stop at a worse point, and nothing in its result would tell you.
The timings show what each route costs. The grey box is a single local NLP; the twin pays for spatial branch-and-bound to get the certificate.
Standard pipe sizes: the discrete version#
Pipe comes in standard sizes, so the real decision is discrete. The obvious formulation gives each segment a binary choice among the sizes. It is a MINLP that contains an opaque block, and discopt’s default global solver refuses it. Without a valid relaxation of the block, branch-and-bound cannot prune soundly:
SIZES = np.array([0.080, 0.100, 0.125, 0.150, 0.200, 0.250, 0.300]) # m
def build_discrete():
m = dm.Model("pipeline_discrete")
D = m.continuous("D", shape=(3,), lb=0.08, ub=0.5)
y = m.binary("y", shape=(3, len(SIZES)))
for i in range(3):
m.subject_to(dm.sum(y[i, :]) == 1)
m.subject_to(D[i] == sum(SIZES[k] * y[i, k] for k in range(len(SIZES))))
f = colebrook(D)
dp = [f[i] * G_DP[i] / D[i] ** 5 for i in range(3)]
m.minimize(sum(A_CAP * L[i] * D[i] ** 1.5 + B_POW * Q[i] * dp[i] for i in range(3)))
m.subject_to(sum(dp) <= DP_MAX)
return m, D
m_disc, D_disc = build_discrete()
try:
m_disc.solve()
except ValueError as exc:
print(f"ValueError: {str(exc)[:260]}...")
Rust model repr unavailable — FBBT/presolve disabled: TypeError: Unknown expression type: CustomCall
ValueError: Model contains a dm.custom(...) AD-only user function that is OUTSIDE the sound reduced-space (MCBox) scope, together with integer/binary variables. Global branch-and-bound needs a valid node relaxation, which a non-MCBox-relaxable opaque callable cannot provi...
The error points to solver="mip-nlp", which is outer approximation
[Duran and Grossmann, 1986]. It accepts opaque blocks, but it cannot do much with this
model, and it is worth seeing why. Its first master proposes the smallest pipe
in every segment, which breaks the pressure budget. For that fixed configuration
the NLP solver detects infeasibility only locally. With no relaxation of the
block, discopt will not promote a local detection into a certified “infeasible”,
so the subproblem ends as error. OA then refuses to add a cut it cannot
justify and stops. The objective also cannot be shown convex through the opaque
block, so OA adds no objective cuts either.
res_oa = m_disc.solve(solver="mip-nlp", time_limit=120)
trace = res_oa.mip_nlp_trace
print(f"status : {res_oa.status} objective: {res_oa.objective}")
print(f"termination reason : {trace['termination_reason']}")
print(f"fixed-NLP statuses : {trace['summary']['fixed_nlp_call_status_counts']}")
OA: nonlinear objective is not convex in the optimization sense; disabling master lower-bound updates and skipping objective OA cuts
No valid dual bound was produced for model 'pipeline_discrete': the relaxation layer could not bound this objective (status=unknown). The result cannot be certified globally optimal. A nonlinear term with no envelope is the usual cause; an epigraph reformulation (minimize z subject to f(x) <= z) often gives the relaxation something to bound.
status : unknown objective: None
termination reason : nonrigorous_nlp_failure
fixed-NLP statuses : {'error': 2}
OA returns unknown with no incumbent. That is the right outcome: it has no
answer it can defend, and it reports none rather than a guess.
A route that works: the continuous optimum plus its neighbours#
The grey box is still an ordinary Python function, and that gives a simple,
transparent route. Take the continuous optimum, list every combination of the
standard sizes just below and just above each diameter (\(2^3 = 8\) candidates),
and evaluate cost and budget directly with friction. Calling the external code
yourself needs no solver.
def cost_and_drop(D):
D = np.asarray(D, dtype=float)
dp = friction(D) * G_DP / D**5
return float(np.sum(A_CAP * L * D**1.5 + B_POW * Q * dp)), float(dp.sum())
neighbours = [(SIZES[SIZES <= d].max(), SIZES[SIZES >= d].min()) for d in D_opt]
print("bracketing sizes per segment:", [tuple(map(float, n)) for n in neighbours])
rows = []
for combo in itertools.product(*neighbours):
cost, drop = cost_and_drop(combo)
rows.append((cost, drop, combo))
rows.sort()
print(f"\n{'D (m)':<24}{'cost $/yr':>12}{'drop bar':>10} ok")
for cost, drop, combo in rows:
ok = "yes" if drop <= DP_MAX else "no"
print(f"{str(tuple(map(float, combo))):<24}{cost:>12,.0f}{drop / 1e5:>10.3f} {ok}")
best_round = next(r for r in rows if r[1] <= DP_MAX)
print(f"\nbest feasible neighbour: {tuple(map(float, best_round[2]))} ${best_round[0]:,.2f}/yr")
bracketing sizes per segment: [(0.15, 0.2), (0.125, 0.15), (0.1, 0.125)]
D (m) cost $/yr drop bar ok
(0.15, 0.125, 0.1) 72,283 4.561 no
(0.15, 0.125, 0.125) 73,758 3.370 no
(0.15, 0.15, 0.1) 74,065 3.939 no
(0.15, 0.15, 0.125) 75,540 2.748 yes
(0.2, 0.125, 0.1) 78,725 3.216 no
(0.2, 0.125, 0.125) 80,200 2.025 yes
(0.2, 0.15, 0.1) 80,506 2.594 yes
(0.2, 0.15, 0.125) 81,982 1.403 yes
best feasible neighbour: (0.15, 0.15, 0.125) $75,540.02/yr
Is the best neighbour the best standard-size design overall? Checking the neighbourhood is a heuristic: a coupled budget can make a combination far from the continuous optimum the cheapest feasible one. With seven sizes and three segments there are only \(7^3 = 343\) designs, so we can enumerate them all and check.
feasible = []
for combo in itertools.product(SIZES, repeat=3):
cost, drop = cost_and_drop(combo)
if drop <= DP_MAX:
feasible.append((cost, combo))
best_enum = min(feasible)
print(f"feasible designs : {len(feasible)} of {len(SIZES) ** 3}")
print(f"enumerated optimum : {tuple(map(float, best_enum[1]))} ${best_enum[0]:,.2f}/yr")
print(f"neighbour heuristic: {tuple(map(float, best_round[2]))} ${best_round[0]:,.2f}/yr")
print(
f"continuous optimum : ${res_ext.objective:,.2f}/yr (a lower bound on any standard design, "
f"certified here by the twin)"
)
assert best_enum[1] == best_round[2]
feasible designs : 112 of 343
enumerated optimum : (0.15, 0.15, 0.125) $75,540.02/yr
neighbour heuristic: (0.15, 0.15, 0.125) $75,540.02/yr
continuous optimum : $72,349.67/yr (a lower bound on any standard design, certified here by the twin)
On this instance the neighbourhood contains the discrete optimum. The certified continuous optimum is a valid lower bound on every standard-size design, because the discrete designs are a subset of the continuous box. So the gap below is the most that restricting to standard sizes can cost on this instance, measured against a bound rather than against another heuristic:
print(
f"standard sizes cost {100 * (best_enum[0] / res_sym.objective - 1):.2f}% more than "
f"the certified continuous optimum"
)
standard sizes cost 4.41% more than the certified continuous optimum
Takeaways#
Use
dm.externalwhen the derivatives belong to someone else’s code, such as a property package, a converged flowsheet, or an iterative solver like this one. Implicit differentiation of the converged residual is the standard way to get those derivatives, and here it took about ten lines.Check the derivatives against finite differences before solving. discopt catches a Jacobian of the wrong shape, but it cannot catch one with the right shape and wrong values.
Keep the opaque part as small as possible. Here only \(f\) is external, and the \(D^{-5}\) law, the costs and the budget are algebraic. Everything outside the block is still visible to discopt.
A grey-box answer is local.
status="feasible"withbound=Nonedescribes exactly what was computed. Where the block has an algebraic form, a symbolic twin can certify a particular instance, as it did here.With discrete decisions, discopt’s solvers will not guess. The global solver refuses the model and OA reports
unknown. Because the block is plain Python, you can still evaluate candidate designs directly, and a certified continuous relaxation gives you a bound to judge them against.