quicopt.stochastic

stochastic

quicopt.stochastic — optimization under uncertainty, written in Pyomo.

Part of a model's data is often unknown when the decision has to be made: demand, prices, yields, arrival times. Declare that data as random variables carrying distributions, and the model is solved over a sample of scenarios drawn from them.

Two rules describe the whole surface:

• A Pyomo variable given a distribution by set_distribution is a random variable, not a decision variable. Every use of the variable references the same sample.
• An expression containing a random variable is itself random, and cannot serve as an objective or a constraint until an aggregator reduces it over the scenarios: expectation for the mean, cvar for the tail, prob for a chance constraint.

set_scenarios sets how many scenarios are drawn and from which seed. Both belong to the model, so repeated solves see the same sample.

Example

Order x units at 3 apiece against demand you will only learn later, pay 10 per unit of mismatch, and meet demand in at least 90% of scenarios:

import pyomo.environ as pyo
from quicopt import Client
from quicopt.stochastic import (Normal, set_distribution, set_scenarios,
                                expectation, prob, maximum)

m = pyo.ConcreteModel()
m.x = pyo.Var(bounds=(0, 200))              # decide now
m.demand = pyo.Var()                        # learn later
set_distribution(m, m.demand, Normal(100.0, 15.0))
set_scenarios(m, 512, seed=42)

m.cost = pyo.Objective(
    expr=3 * m.x + 10 * expectation(maximum(m.demand - m.x, 0)),
    sense=pyo.minimize)
m.service = pyo.Constraint(expr=prob(m.demand - m.x <= 0) >= 0.9)

result = Client().solve(m)
print(result.solution)                      # {'x1': 107.9…}

Every mention of m.demand refers to the same sample. Two independent random variables are two Pyomo variables under two names.

maximum and minimum build the kinked expressions that price recourse: maximum(demand - stock, 0) is the shortfall, and zero when there is none.

Models written this way are solved by the service. Attaching a local solver to one fails — an expectation has no value until the scenarios are drawn.

Empirical dataclass

Empirical(data: list)

A random variable given as a fixed scenario column of data.

Exactly Program.scenarios values, one per scenario. Several empirical columns are read at the same scenario index, so columns drawn jointly stay correlated — which is how a joint distribution is expressed.

Distribution

Distribution(head, *params)

Any distribution the service supports, named by head.

Use :class:Normal where it fits; use this for a distribution that has no named class here yet, as in Distribution("uniform", 0.0, 1.0).

A parameter may be a number or a deterministic Pyomo expression. An expression gives an endogenous distribution — one whose parameters depend on the decision, such as a demand whose mean rises with the price you set.

Draw from the distribution head with the given parameters, each a number or a deterministic Pyomo expression.

Source code in src/quicopt/stochastic.py

def __init__(self, head, *params):
    """Draw from the distribution ``head`` with the given parameters, each a
    number or a deterministic Pyomo expression."""
    self.head = str(head)
    self.params = params

Normal

Normal(mu, sigma)

Bases: Distribution

The normal distribution with mean mu and standard deviation sigma.

Declare a draw from N(mu, sigma) — the mean and the standard deviation, each a number or a deterministic Pyomo expression.

Source code in src/quicopt/stochastic.py

def __init__(self, mu, sigma):
    """Declare a draw from ``N(mu, sigma)`` — the mean and the standard
    deviation, each a number or a deterministic Pyomo expression."""
    super().__init__("normal", mu, sigma)

set_distribution

set_distribution(model, v, dist, name=None)

Turn the Pyomo variable v into a random variable drawn from dist.

v is no longer a decision variable: the solver is handed its value rather than choosing it, and every use of v means the same sample within a scenario.

dist is a :class:Distribution such as Normal(100.0, 15.0), or an :class:~quicopt.ir.Empirical column holding one observed value per scenario. name defaults to the variable's own name and is what identifies the random variable, so two independent ones need two names. Returns v.

Give v no bounds and no domain. Its distribution already says what values it takes, and a variable carrying both is rejected on import.

Raises:

Type	Description
TypeError	If v is not a Pyomo variable, or dist is neither a Distribution nor an Empirical column.
ValueError	If name is already taken by a different random variable.

Source code in src/quicopt/stochastic.py

def set_distribution(model, v, dist, name=None):
    """Turn the Pyomo variable ``v`` into a random variable drawn from ``dist``.

    ``v`` is no longer a decision variable: the solver is handed its value rather
    than choosing it, and every use of ``v`` means the same sample within a scenario.

    ``dist`` is a :class:`Distribution` such as ``Normal(100.0, 15.0)``, or an
    :class:`~quicopt.ir.Empirical` column holding one observed value per scenario.
    ``name`` defaults to the variable's own name and is what identifies the random
    variable, so two independent ones need two names. Returns ``v``.

    Give ``v`` no bounds and no domain. Its distribution already says what values
    it takes, and a variable carrying both is rejected on import.

    Raises:
        TypeError: If ``v`` is not a Pyomo variable, or ``dist`` is neither a
            ``Distribution`` nor an ``Empirical`` column.
        ValueError: If ``name`` is already taken by a different random variable.
    """
    if not (hasattr(v, "is_variable_type") and v.is_variable_type()):
        raise TypeError(f"a distribution attaches to a Pyomo variable, not to {type(v).__name__}")
    if not isinstance(dist, (Distribution, Empirical)):
        raise TypeError("a distribution must be a Distribution (e.g. Normal(mu, sigma)) "
                        f"or an Empirical column, got {type(dist).__name__}")

    record = _declare(model)
    declared = str(v.name if name is None else name)
    for (other, other_name, _) in record.sources.values():
        if other is not v and other_name == declared:
            raise ValueError(f"the name '{declared}' already belongs to another random variable — "
                             "one name is one random variable, so independent ones need "
                             "distinct names")
    record.sources[id(v)] = (v, declared, dist)
    return v

set_scenarios

set_scenarios(model, n, seed=None)

Draw n scenarios, optionally from a given seed.

More scenarios estimate the true problem more closely and cost more to solve. Both settings belong to the model, not to the solve, so the same model always faces the same sample and two solves of it are comparable. Left unset, a model is solved over one scenario, where every random variable takes a single value.

n and seed are both at least 1. Leaving seed out keeps the current one.

Raises:

Type	Description
ValueError	If n or seed is below 1.

Source code in src/quicopt/stochastic.py

def set_scenarios(model, n, seed=None):
    """Draw ``n`` scenarios, optionally from a given ``seed``.

    More scenarios estimate the true problem more closely and cost more to solve.
    Both settings belong to the model, not to the solve, so the same model always
    faces the same sample and two solves of it are comparable. Left unset, a model
    is solved over one scenario, where every random variable takes a single value.

    ``n`` and ``seed`` are both at least 1. Leaving ``seed`` out keeps the current
    one.

    Raises:
        ValueError: If ``n`` or ``seed`` is below 1.
    """
    # 0 is rejected rather than defaulted: protobuf cannot tell a zero from an absent
    # field, so a 0 sent here would arrive at the service as "use your default".
    if n < 1:
        raise ValueError(f"a model is solved over at least 1 scenario, got {n}")
    record = _declare(model)
    record.scenarios = int(n)
    if seed is not None:
        if seed < 1:
            raise ValueError(f"the scenario seed must be at least 1, got {seed}")
        record.seed = int(seed)

expectation

expectation(x)

The expected value of x over the scenarios.

Minimizing an expectation optimizes the average case and says nothing about the bad ones; use :func:cvar when the bad ones are what matter.

x is any expression containing a random variable. The result is deterministic, and can be used anywhere a number can.

Source code in src/quicopt/stochastic.py

def expectation(x):
    """The expected value of ``x`` over the scenarios.

    Minimizing an expectation optimizes the average case and says nothing about
    the bad ones; use :func:`cvar` when the bad ones are what matter.

    ``x`` is any expression containing a random variable. The result is
    deterministic, and can be used anywhere a number can.
    """
    return _Aggregate("smean", [x])

cvar

cvar(x, alpha)

The conditional value at risk of x at level alpha.

The mean of x over its worst 1 − alpha fraction of scenarios — at alpha = 0.95, the average of the worst 5%. Minimizing it optimizes the tail instead of the average, and is the usual way to ask for a solution that holds up in bad scenarios rather than merely on average.

x is any expression containing a random variable. alpha is the tail level, a plain number strictly between 0 and 1; it cannot depend on a decision.

Raises:

Type	Description
TypeError	If alpha is not a plain number.
ValueError	If alpha does not lie strictly between 0 and 1.

Source code in src/quicopt/stochastic.py

def cvar(x, alpha):
    """The conditional value at risk of ``x`` at level ``alpha``.

    The mean of ``x`` over its worst ``1 − alpha`` fraction of scenarios — at
    ``alpha = 0.95``, the average of the worst 5%. Minimizing it optimizes the
    tail instead of the average, and is the usual way to ask for a solution that
    holds up in bad scenarios rather than merely on average.

    ``x`` is any expression containing a random variable. ``alpha`` is the tail
    level, a plain number strictly between 0 and 1; it cannot depend on a decision.

    Raises:
        TypeError: If ``alpha`` is not a plain number.
        ValueError: If ``alpha`` does not lie strictly between 0 and 1.
    """
    if not _is_number(alpha):
        raise TypeError(f"the tail level must be a plain number, got {type(alpha).__name__}")
    if not 0.0 < alpha < 1.0:
        raise ValueError(f"the tail level must lie strictly between 0 and 1, got {alpha}")
    return _Aggregate("scvar", [x, float(alpha)])

prob

prob(relation)

The probability that relation holds — the fraction of scenarios it does.

This is what a chance constraint is built from::

m.service = pyo.Constraint(expr=prob(m.demand - m.x <= 0) >= 0.9)

which reads as demand is met in at least 90% of scenarios. Note that the line holds two comparisons: the one inside prob is the event being measured, the outer one is the service level demanded of it.

relation is a comparison, a <= b or a >= b, with at least one side containing a random variable. The result is a probability between 0 and 1.

Raises:

Type	Description
TypeError	If the argument is not a comparison. An equality is refused too: for a continuous quantity its probability is zero.
ValueError	If the comparison is strict. Use <= or >=, which for a continuous quantity mean the same thing anyway.

Source code in src/quicopt/stochastic.py

def prob(relation):
    """The probability that ``relation`` holds — the fraction of scenarios it does.

    This is what a chance constraint is built from::

        m.service = pyo.Constraint(expr=prob(m.demand - m.x <= 0) >= 0.9)

    which reads as *demand is met in at least 90% of scenarios*. Note that the
    line holds two comparisons: the one inside ``prob`` is the event being
    measured, the outer one is the service level demanded of it.

    ``relation`` is a comparison, ``a <= b`` or ``a >= b``, with at least one side
    containing a random variable. The result is a probability between 0 and 1.

    Raises:
        TypeError: If the argument is not a comparison. An equality is refused
            too: for a continuous quantity its probability is zero.
        ValueError: If the comparison is strict. Use ``<=`` or ``>=``, which for a
            continuous quantity mean the same thing anyway.
    """
    if not isinstance(relation, InequalityExpression):
        raise TypeError("prob takes a comparison, as in prob(demand - x <= 0), "
                        f"not {type(relation).__name__}")
    if relation.strict:
        raise ValueError("prob takes <= or >=; a strict comparison has no scenario "
                         "counterpart (and for a continuous quantity means the same thing)")

    lower, upper = relation.args              # Pyomo normalizes any inequality to lower <= upper
    if _is_number(upper):
        return _Aggregate("sfreq_leq", [lower, float(upper)])
    if _is_number(lower):
        return _Aggregate("sfreq_geq", [upper, float(lower)])
    return _Aggregate("sfreq_leq", [lower - upper, 0.0])   # P(a ≤ b) = P(a − b ≤ 0)

maximum

maximum(*args)

The largest of two or more expressions.

This is how recourse is priced: maximum(demand - stock, 0) is the shortfall, and zero whenever stock covers demand. Python's built-in max cannot do this — it would have to compare the arguments, and their values are not known yet.

Raises:

Type	Description
ValueError	If given fewer than two arguments.

Source code in src/quicopt/stochastic.py

def maximum(*args):
    """The largest of two or more expressions.

    This is how recourse is priced: ``maximum(demand - stock, 0)`` is the
    shortfall, and zero whenever stock covers demand. Python's built-in ``max``
    cannot do this — it would have to compare the arguments, and their values are
    not known yet.

    Raises:
        ValueError: If given fewer than two arguments.
    """
    return _ne.MaxExpression(_at_least_two("maximum", args))

minimum

minimum(*args)

The smallest of two or more expressions.

The counterpart of :func:maximum — minimum(demand, stock) is what you actually sell.

Raises:

Type	Description
ValueError	If given fewer than two arguments.

Source code in src/quicopt/stochastic.py

def minimum(*args):
    """The smallest of two or more expressions.

    The counterpart of :func:`maximum` — ``minimum(demand, stock)`` is what you
    actually sell.

    Raises:
        ValueError: If given fewer than two arguments.
    """
    return _ne.MinExpression(_at_least_two("minimum", args))