Uncertainty#

PULPO optimizes one impact under the uncertainty its data declare. The method needs no sampling to solve: the impact’s mean and variance follow in closed form, the chance constraints become a second-order cone, and the cone is solved in reduced space, over the choice alternatives only. Sampling enters afterwards, to validate a solution out of sample. The uncertainty notebook runs every step on the bundled demo database with open-source solvers.

Note

As of now, uncertainty in the LCA data is considered only in the biosphere (intervention) flows and the characterization factors. Technosphere exchanges are deterministic, even where the database declares a distribution for them (assumption A4 below). Uncertain capacities or availabilities enter as uncertain process bounds.

from pulpo.utils import uncertainty as unc

data  = unc.import_declared(worker)                         # parameters of the impact
unc.override(data, 'If', 'my_foreground_db', expert_specs)  # elicited distributions
mom   = unc.compute_moments(data, worker)                   # closed-form mean and variance
ccp   = unc.ChanceConstrained(worker, mom, upper_bounds={activity: capacity_spec})
front = ccp.solve([0.5, 0.9, 0.99])                         # reliability sweep
dec   = unc.decompose(front[0.9], mom)                      # exact Sobol' indices
scr   = unc.screen_undeclared(front[0.9], data, worker, exact_cfs=unc.co2_flows(worker))
val   = unc.validate(front, ccp, data, n=200_000, seed=1)

pulpo.pulpo_unc.PulpoOptimizerUnc keeps data on the worker and offers the import, the override, the moments, the chance-constrained problem, the screening and the validation as methods (import_uncertainty_data, apply_expert_knowledge, moments, chance_constrained, screen_undeclared, validate). decompose and the other analyses are called as above, with worker.uncertainty_data as the data.

Data and assumptions#

\(A\) is the technosphere matrix (square, invertible, deterministic), \(B\) the biosphere matrix with entries \(b_{ej}\), \(q\) the characterization factors (CFs) of one LCIA method, \(s\) the scaling vector. The impact is \(X(s) = \sum_e q_e \sum_j b_{ej} s_j\).
The uncertain inputs are the entries of \(B\) and the CFs that declare a distribution (in the database, in the method, or by expert override), and the bounds on processes that the problem declares uncertain (e.g. an availability or a capacity).
Supported families: lognormal (also for negative amounts), normal, uniform, triangular, and exact values (uncertainty_type 1). Any other family raises.

The method rests on five assumptions, each stated where it is used:

	Assumption	Where
A1	all uncertain inputs are mutually independent	moments, decomposition
A2	\(X(s)\) is approximately normal (the only normality assumption)	the impact’s chance constraint; tested by the validation
A3	undeclared parameters are deterministic	everywhere; tested by the screening
A4	\(A\) is deterministic	makes the reduced space exact
A5	the failure budget is split by Boole’s inequality	joint chance constraints; conservative under any dependence

import_declared takes every nonzero entry of \(B\) on a characterized flow and every CF on a flow that occurs in \(B\). Parameters without a distribution are kept as undeclared: their mean is their amount and their variance zero. No contribution filter is applied and nothing is gap-filled.

Moments#

\[\begin{split} \begin{aligned} \mathrm{E}[X] &= \mu^\top s, & \mu_j &= \textstyle\sum_e \mathrm{E}[q_e]\,\mathrm{E}[b_{ej}] \\ \mathrm{Var}\,X &= \textstyle\sum_j d_j s_j^2 + \sum_e w_e y_e^2, & d_j &= \textstyle\sum_e (\mathrm{E}[q_e]^2 + \mathrm{Var}\,q_e)\,\mathrm{Var}\,b_{ej},\quad w_e = \mathrm{Var}\,q_e \end{aligned} \end{split}\]

with \(y_e = \sum_j \mathrm{E}[b_{ej}] s_j\) the mean flow \(e\). The second sum is the covariance that every process emitting flow \(e\) shares through its one factor \(q_e\); Moments.std_independent drops it for comparison.

Chance constraints#

At a joint reliability level \(\lambda\) over \(K\) events (the impact row and one per uncertain bound), each event may fail with probability \(\varepsilon_k = w_k(1-\lambda)\), equal weights by default. Boole’s inequality then bounds the probability that any event fails by \(1-\lambda\):

\[\begin{split} \begin{aligned} \text{impact row:}\quad & P(X \le z) \ge \lambda_z = 1-\varepsilon_0 &&\Rightarrow\quad \mu^\top s + \kappa\,\sigma(s) \le z,\quad \kappa = \Phi^{-1}(\lambda_z) \\ \text{bound } j:\quad & P(s_j \le U_j) \ge 1-\varepsilon_k &&\Rightarrow\quad s_j \le F_{U_j}^{-1}(\varepsilon_k) \end{aligned} \end{split}\]

and minimizing \(z\) gives

\[ \min_s\ \mu^\top s + \kappa\,\sigma(s)\quad \text{s.t. every constraint of the instance, with the uncertain bounds at their quantiles.} \]

The bounds enter at the exact quantile of their declared family; only the impact row uses the normal approximation (A2). \(\kappa \ge 0\) (\(\lambda_z \ge 1/2\)) keeps the problem convex; lower levels are refused. allocation='individual' imposes every event at \(\lambda\) on its own, for comparison: it controls no joint probability.

Reduced space#

With \(A\) square and invertible, every scaling vector that meets the balances is

\[ s = s_0 + S v,\qquad s_0 = A^{-1}\tilde f,\qquad S = A^{-1}E, \]

with one free variable \(v_k\) per alternative (the net output on its product row) and \(\tilde f\) the demand on all other rows. Every static PULPO constraint is linear in \(s\), so it becomes one row in \(v\), and the problem over \(v\) is the problem over \(s\): the reduction is exact. \(S\) is never formed; each row \(m^\top S\) costs one adjoint solve with one factorization of \(A\) (PARDISO, UMFPACK or SciPy’s SuperLU). Every finite process bound is a row \(S[j, :]\), so the problem stays small as long as only the processes with a real limit are bounded: finite default_limits on every process would put the whole of \(S\) into it.

The deterministic LP is solved this way by solve(formulation='reduced'). For the chance-constrained problem \(\sigma(s) = \lVert R\,[1; v]\rVert\) with \(R^\top R = G^\top G\), \(G = Q^{1/2}[s_0, S]\) and \(Q = \mathrm{diag}(d) + B_u^\top \mathrm{diag}(w) B_u\): a second-order cone with one column per alternative, solved by Clarabel or Gurobi.

Exact variance decomposition#

\(X\) is a sum of products of independent inputs, so its Sobol’ decomposition stops at second order. At a fixed \(s\), with \(V = \mathrm{Var}\,X\):

\[\begin{split} \begin{aligned} S_1(q_e) &= w_e y_e^2 / V, & S_1(b_{ej}) &= \mathrm{E}[q_e]^2\,\mathrm{Var}(b_{ej})\,s_j^2 / V, & S_2(q_e, b_{ej}) &= w_e\,\mathrm{Var}(b_{ej})\,s_j^2 / V, \\ S_T(q_e) &= w_e \big(y_e^2 + \textstyle\sum_j \mathrm{Var}(b_{ej})\,s_j^2\big) / V, & S_T(b_{ej}) &= (\mathrm{E}[q_e]^2 + w_e)\,\mathrm{Var}(b_{ej})\,s_j^2 / V \end{aligned} \end{split}\]

and \(\sum S_1 + \sum S_2 = 1\) exactly. decompose reports them per parameter and per family (e.g. per database with families_by_database, and the CFs). B entries never interact with each other, nor CFs with each other, so a family’s indices are the sums over its members. The uncertain bounds are not inputs of \(X\) at a fixed \(s\); they shape \(s\).

Undeclared parameters#

screen_undeclared tests assumption A3. Every undeclared parameter gets the same coefficient of variation \(r\) with its mean unchanged (\(\mathrm{Var}\,b = r^2 b_0^2\), \(\mathrm{Var}\,q = r^2\,\mathrm{E}[q]^2\)) and the indices above are evaluated at a fixed decision. The CFs passed as exact_cfs never receive a width: for a global warming potential these are the CO₂ flows, whose factor is exact by definition because CO₂ is the reference gas (co2_flows finds them by name in an ecoinvent biosphere). Every undeclared B entry is widened: none is exact by definition, holding one exact could only lower the figures, and only the few at the top of the ranking need an expert’s judgement. The parameters are ranked by their deterministic contribution \(\mathrm{E}[q_e]\,b_{0,ej}\,s_j\), which fixes the order of their total-order indices up to the CF factor \(1 + w_e/\mathrm{E}[q_e]^2\) and does not depend on \(r\); the indices at each \(r\) are reported beside it.

width_sensitivity evaluates \(\sigma(s^*)\) with separate widths per database and for the CFs. The mean does not change, so \(\kappa\,\Delta\sigma\) is an upper bound on the rise of the optimal adjusted impact (re-optimizing can only do better). widen(data, r, exact_cfs=co2_flows(worker)) returns the widened data, with a lognormal of the same mean and coefficient of variation (mirrored for negative amounts), so the same data give a re-solved front and validation draws. Widths are given per database name and 'Cf'; a name that is not a subgroup of the data raises.

Out-of-sample validation#

validate holds each solution \(s^*\) and its target \(z^*\) fixed and draws every declared input from its own family, together with the uncertain bounds, from one recorded seed. Every decision is priced on the same draws, so comparisons between decisions are paired. Per point it reports:

the coverage of the impact row against \(\lambda_z\), of each bound against \(1-\varepsilon_k\), and of all events together against \(\lambda\), with Wilson intervals;
the normal-approximation error (the empirical \(\lambda_z\) quantile minus \(z^*\)), the skewness and excess kurtosis of \(X\);
the empirical CVaR at \(\lambda_z\) against the Gaussian \(\mu + \varphi(\Phi^{-1}(\lambda_z))/(1-\lambda_z)\,\sigma\);
the regret against the lowest level of the front and against the best front point in each draw (mean and 95th percentile).

Further decisions, such as the deterministic optimum, are priced against the same targets with designs=. Each parameter is drawn from its own stream, keyed by the seed and its position, so its draws depend on the seed alone; decisions validated in separate calls with the same seed are paired as well. Only B entries that carry variance at the decisions are drawn; the smallest ones, below tol of each decision’s variance together, are held at their means, and the share drawn is reported.

Diagnostics#

diagnostics(front, moments) reports at every point: \(\sigma\) without the shared-CF covariance and the ratio \(\sigma/\sigma_{\text{indep}}\); the CFs’ share of the variance; the bound imposed on each uncertain process; the size of the reduced problem and the solve times.

Building blocks#

The reduced system and the sampler are public, so a study can formulate its own problem over \(v\), such as a scenario-based CVaR front:

reduced.build(worker) returns the ReducedModel of an instance: its ReducedSystem maps \(v\) to \(s\) (base, recover) and projects any linear functional of \(s\) (project, rows, project_vectors); linear_program() assembles the constraint rows in \(v\).
ChanceConstrained.projections() returns \(s_0\), \(\mathrm{E}[X] = m_0 + m^\top v\), the rows \(S_J\) of \(S\) on the processes that carry variance and \(\mathrm{E}[B_u]S\).
draw_parameters returns draws of the declared parameters as arrays, and sample_specs draws any list of specs.

Solvers and licences#

Problem	Solvers	Licence
deterministic LP, `formulation='full'`	HiGHS (default), Gurobi, GAMS/NEOS	HiGHS: none
deterministic LP, `formulation='reduced'`	HiGHS (default), Gurobi	HiGHS: none
chance-constrained cone	Clarabel (default), Gurobi	Clarabel: none

HiGHS and Clarabel are installed with PULPO and need no licence. Gurobi is used when gurobipy is installed; the size-limited licence that comes with pip install gurobipy is for non-production use and covers reduced problems of up to 2,000 variables and 2,000 linear constraints, or 200 variables once quadratic terms (the cone) are present (see Gurobi’s licence terms). A reduced problem has one column per alternative, so most case studies fit.

The defaults stay formulation='full' and instantiate(scale=False). Current development considers formulation='reduced' and scale=True superior, and a future release may switch the defaults.