Contents

Example finding multiple solutions

Example finding multiple solutions

This example shows how multiple solutions can be found when the model is under-constrained.

Load example model

A PACTOT material flow model is defined in example_model.py. This was chosen to be fairly simple, but complicated enough to demonstrate enough features of the solver. The system describes the steelmaking, construction and waste processing.

Load and solve the model, solution vector X and flows table:

import numpy as np
import pandas as pd

from probs_models.observations import (
    ProductionObs,
    ImportsObs,
    ExportsObs,
    ProcessInputObs,
)
from probs_models.models import pactot

import example_model

def solve(observations):
    m = pactot.build_model(example_model.recipe_builders, 
                           example_model.objects,
                           observations)
    pactot.diagnose(m)
    result = pactot.solve(m, max_solutions=None)
    return result

def get_flows(result):
    flows = [pactot.build_flows_table(result.model, X) for X in result.solutions]
    for f, label in zip(flows, result.solution_labels):
        f.loc[f["value"] < 5e-2, "value"] = 0
        f.loc[:, "solution"] = label
    return pd.concat(flows, axis=0, ignore_index=True)

result = solve([
    ProductionObs(["Buildings"], value=10),
    ImportsObs(["Steel"], value=1, quality=1),
    ExportsObs(["Steel"], value=0, quality=1),
    ExportsObs(["Coal"], value=0, quality=1),
    ImportsObs(["CO2"], value=0, quality=1),
    ImportsObs(["Scrap steel"], value=0, quality=1),
    ImportsObs(["Scrap wood"], value=0, quality=1),
    ImportsObs(["Mixed waste"], value=0, quality=1),
])
flows = get_flows(result)
print(result.status, result.solution_labels)

SolutionStatus.UNIQUE ['First solution']

Show the basic Sankey diagram (using the layout defined in sdd.py):

from sdd import sdd
from floweaver import weave, QuantitativeScale, CategoricalScale

weave(sdd, flows).to_widget(debugging=True, width=1000, height=300)

Under-constrained solution

If we repeat without any observation of exports of Steel, the solution is under-constrained. In this case 4 example alternative solutions have been returned; the labels explained what was being adjusted in each case.

result = solve([
    ProductionObs(["Buildings"], value=10),
    ImportsObs(["Steel"], value=1, quality=1),
    #ExportsObs(["Steel"], value=0, quality=1),
    ExportsObs(["Coal"], value=0, quality=1),
    ImportsObs(["CO2"], value=0, quality=1),
    ImportsObs(["Scrap steel"], value=0, quality=1),
    ImportsObs(["Scrap wood"], value=0, quality=1),
    ImportsObs(["Mixed waste"], value=0, quality=1),
])
flows = get_flows(result)
print(result.status, result.solution_labels)

SolutionStatus.UNDERCONSTRAINED ['Push s_ProducingIronOre x 0', 'Push s_ProducingOxygen x 2', 'Push s_Ironmaking x 2', 'Push w_Steel x 2']

For example, see the possible solutions for the flows related to Steel:

flows.query("material == 'Steel'") \
    .set_index(["source", "target", "solution"])["value"] \
    .unstack()
solution Push s_Ironmaking x 2 Push s_ProducingIronOre x 0 Push s_ProducingOxygen x 2 Push w_Steel x 2
source target
Steel Fabrication 6.250000 6.25 6.250 6.250000
exports Steel 6.994575 0.00 7.875 1.744575
Steelmaking Steel 12.244575 5.25 13.125 6.994575
imports Steel Steel 1.000000 1.00 1.000 1.000000

It’s easier to see the solution graphically:

def show_sankey_solution_variation(flows):
    grouped = (
        flows.groupby(["source", "target", "material"], as_index=False)["value"]
        .agg({"value": "mean", "std": "std"})
    )
    return (
        weave(sdd, 
              grouped,
              link_color=QuantitativeScale("std", intensity="value", domain=(0,1)),
              link_width="value",
              measures=["value", "std"])
        .to_widget(debugging=True, width=1000, height=300)
    )

show_sankey_solution_variation(flows)

/Users/rcl38/work/probs-models/.venv/lib/python3.9/site-packages/floweaver/color_scales.py:122: RuntimeWarning: invalid value encountered in double_scalars
  value /= measures[self.intensity]

In this diagram the colour intensity shows the variation between the alternative solutions (dark red means that part of the system is under-constrained).

Since the actual mean value returned is fairly arbitrary, maybe we should just pick one of the alternative solutions:

def show_one_sankey_solution(flows):
    scale = flows.groupby("solution")["value"].sum()
    smallest_solution = scale.index[scale.argmin()]
    variation = (
        flows.groupby(["source", "target", "material"])
        ["value"].std()
    )
    flows = flows[flows["solution"] == smallest_solution].set_index(["source", "target", "material"])
    flows["var"] = variation > 1e-3
    return (
        weave(sdd, 
              flows.reset_index(),
              link_color=CategoricalScale("var", {False: "lightsteelblue", True: "tomato"}),
              link_width="value",
              measures=["value", "var"])
        .to_widget(debugging=True, width=1000, height=300)
    )

show_one_sankey_solution(flows)

Another example under-constrained solution

Another example:

result = solve([
    ProductionObs(["Buildings"], value=10),
    # New:
    ProductionObs(["Iron"], value=3),
    #ImportsObs(["Steel"], value=1, quality=1),
    #ExportsObs(["Steel"], value=0, quality=1),
    ExportsObs(["Coal"], value=0, quality=1),
    ImportsObs(["CO2"], value=0, quality=1),
    ImportsObs(["Scrap steel"], value=0, quality=1),
    ImportsObs(["Scrap wood"], value=0, quality=1),
    ImportsObs(["Mixed waste"], value=0, quality=1),
])
flows = get_flows(result)
result.solution_labels

['Push v_Steel x 0', 'Push w_Steel x 2']

show_one_sankey_solution(flows)

Here we have only the trade flows in Steel being under-constrained.

Both under- and over-constrained solution

Finally, let’s demonstrate that the solution can be over-constrained (with observations needing to be reconciled) in places, and under-constrained in others.

result = solve([
    ProductionObs(["Buildings"], value=10),
    ProductionObs(["Iron"], value=3),
    ProductionObs(["Iron ore"], value=3),
    #ImportsObs(["Steel"], value=1, quality=1),
    #ExportsObs(["Steel"], value=0, quality=1),
    ExportsObs(["Coal"], value=0, quality=1),
    ImportsObs(["CO2"], value=0, quality=1),
    ImportsObs(["Scrap steel"], value=0, quality=1),
    ImportsObs(["Scrap wood"], value=0, quality=1),
    ImportsObs(["Mixed waste"], value=0, quality=1),
])
flows = get_flows(result)
result.solution_labels

/Users/rcl38/work/probs-models/.venv/lib/python3.9/site-packages/cvxpy/problems/problem.py:1333: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  warnings.warn(

['Push v_Steel x 0', 'Push w_Steel x 2']

show_one_sankey_solution(flows)

previous

PACTOT model

next

SCOUT content flow model