Skip to content

Optimization API

Used in the optimization guide, resource-limited optimization, and supply-chain examples.

Scenario dataclass 

Scenario(
    objective: str = "automation",
    min_goods_balance: float | None = 0.0,
    banned_pms: tuple[str, ...] = (),
    building_level_caps: tuple[tuple[str, float], ...] = (),
    banned_building_groups: tuple[str, ...] = (),
    throughput_multipliers: tuple[
        tuple[str, float], ...
    ] = (),
    era_cap: int | None = None,
    construction_cost_cap: float | None = None,
    employment_cap: float | None = None,
    arable_land_cap: float | None = None,
    min_net_infrastructure: float | None = None,
    urbanization_per_center: float | None = None,
    name: str | None = None,
    goods_import: tuple[tuple[str, float], ...] = (),
    goods_export: tuple[tuple[str, float], ...] = (),
    workforce_ratio: float = 1.0,
)

An immutable optimization recipe expressed as data.

The optimizer interpreting a scenario owns all economy-dependent compilation. This object records assumptions interpreted by the nominal, market, and equilibrium formulations.

min_goods_balance sets the minimum weekly balance for every good: building output plus imports minus building inputs and exports. Household consumption is excluded. The default of 0.0 requires nonnegative balances; negative floors allow deficits, and None disables the floor.

min_net_infrastructure sets minimum infrastructure generation minus consumption, weighted by building levels. 0.0 prevents a deficit; positive minima require spare capacity, negative minima allow a bounded deficit, and None leaves infrastructure unconstrained.

goods_import and goods_export record fixed weekly (good, amount) pairs with non-negative amounts. Imports add sell orders; exports add buy orders. Repeated entries are summed, and a good may be both imported and exported. These orders enter the goods-balance constraint when enabled.

workforce_ratio is the uniform worker share of total population. It defaults to 1.0; employment constraints still count only workers.

Used by:

Methods:

  • display_name –

    Return the explicit name or one derived from the export goods basket.

display_name 

display_name() -> str

Return the explicit name or one derived from the export goods basket.

Source code in src/vic3_analysis/optimize/scenario.py 
72
73
74
75
76
77
78
def display_name(self) -> str:
    """Return the explicit name or one derived from the export goods basket."""
    if self.name is not None:
        return self.name
    if self.goods_export:
        return "+".join(dict.fromkeys(good for good, _amount in self.goods_export))
    return "unnamed"

BaseOptimizer 

BaseOptimizer(model: Economy)

Bases: ABC

Abstract base for economy-bound optimization formulations.

Subclassed by:

Methods:

  • compile –

    Compile scenario into the concrete optimizer's problem type.

  • summary –

    Return a fresh one-row report of the latest returned solution.

Attributes:

  • goods_input_matrix (ndarray) –

    Return the current problem's adjusted input matrix.

  • goods_output_matrix (ndarray) –

    Return the current problem's adjusted output matrix.

  • goods_matrix (ndarray) –

    Return the current problem's adjusted net-goods matrix.

  • gdp_vector (ndarray) –

    Return the current problem's nominal GDP vector.

  • throughput_multipliers (ndarray) –

    Return the current problem's per-configuration multipliers.

Source code in src/vic3_analysis/optimize/base.py 
163
164
165
166
167
168
169
def __init__(self, model: Economy) -> None:
    """Bind the optimizer to *model*."""
    self.model = model
    self.result = None
    self.scenario = None
    self.problem = None
    self._state = None

goods_input_matrix property 

goods_input_matrix: ndarray

Return the current problem's adjusted input matrix.

goods_output_matrix property 

goods_output_matrix: ndarray

Return the current problem's adjusted output matrix.

goods_matrix property 

goods_matrix: ndarray

Return the current problem's adjusted net-goods matrix.

gdp_vector property 

gdp_vector: ndarray

Return the current problem's nominal GDP vector.

throughput_multipliers property 

throughput_multipliers: ndarray

Return the current problem's per-configuration multipliers.

compile abstractmethod 

compile(scenario: Scenario) -> ProblemT

Compile scenario into the concrete optimizer's problem type.

Source code in src/vic3_analysis/optimize/base.py 
171
172
173
@abstractmethod
def compile(self, scenario: Scenario) -> ProblemT:
    """Compile *scenario* into the concrete optimizer's problem type."""

summary 

summary() -> DataFrame

Return a fresh one-row report of the latest returned solution.

GDP is reported weekly and annually; gdp_per_capita is annual. Employment counts workers, population includes dependents, and wealth is a population-weighted tier. Infrastructure usage is consumption minus generation. Construction cost and arable-land consumption use the model's existing accounting. throughput_multiplier is NaN unless all compiled configuration multipliers are equal.

Reporting never solves or evaluates the problem again. Binding a new problem clears the retained solution. Equilibrium searches that return a verified state on budget exhaustion remain reportable even when search_success is false.

Raises:

  • ValueError –

    If the current problem has no returned solution.

Source code in src/vic3_analysis/optimize/base.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
def summary(self) -> pd.DataFrame:
    """Return a fresh one-row report of the latest returned solution.

    GDP is reported weekly and annually; ``gdp_per_capita`` is annual.
    Employment counts workers, population includes dependents, and wealth
    is a population-weighted tier. Infrastructure usage is consumption
    minus generation. Construction cost and arable-land consumption use
    the model's existing accounting. ``throughput_multiplier`` is NaN
    unless all compiled configuration multipliers are equal.

    Reporting never solves or evaluates the problem again. Binding a new
    problem clears the retained solution. Equilibrium searches that return
    a verified state on budget exhaustion remain reportable even when
    ``search_success`` is false.

    Raises:
        ValueError: If the current problem has no returned solution.
    """
    state = self._require_state()
    problem = self._require_problem()
    result = self.result
    if result is None:
        raise ValueError("Optimizer has no returned solution to summarize.")
    scenario = problem.scenario
    multipliers = problem.throughput_multipliers
    values: dict[str, object] = {
        "objective": scenario.objective,
        "employment_cap": scenario.employment_cap,
        "workforce_ratio": scenario.workforce_ratio,
        "banned_building_groups": ";".join(scenario.banned_building_groups),
        "min_net_infrastructure": scenario.min_net_infrastructure,
        "min_goods_balance": scenario.min_goods_balance,
        "configuration_count": len(self.model.building_index),
        "throughput_multiplier": (
            float(multipliers[0])
            if multipliers.size and np.all(multipliers == multipliers[0])
            else float("nan")
        ),
        "economy_of_scale_level_cap": 0.0,
        "gdp_weekly": state.gdp_weekly,
        "annual_gdp": state.gdp(annual=True),
        "employment": state.total_employment,
        "population": state.total_population,
        "gdp_per_capita": state.gdp_per_capita(annual=True),
        "average_wealth": state.average_wealth,
        "construction_cost": self.model.construction_cost(state),
        "arable_land_consumption": self.model.arable_land_consumption(state),
        "net_infrastructure_usage": float(
            state.building_levels @ self.model.infrastructure_vector
        ),
        "search_success": bool(result.success),
        "search_status": int(result.status),
        "search_message": str(result.message),
    }
    if "nfev" in result:
        values["trial_evaluations"] = int(result.nfev)
    return pd.DataFrame([values])

LinearProblem dataclass 

LinearProblem(
    economy: Economy,
    scenario: Scenario,
    throughput_multipliers: ndarray,
    input_matrix: ndarray,
    output_matrix: ndarray,
    net_matrix: ndarray,
    gdp_vector: ndarray,
    goods_import: ndarray,
    goods_export: ndarray,
    A_ub: ndarray | None,
    b_ub: ndarray | None,
    A_eq: ndarray | None,
    b_eq: ndarray | None,
    bounds: Bounds,
    inequality_slices: ConstraintSlices,
    equality_slices: ConstraintSlices,
    objective_vector: ndarray,
)

Bases: _CompiledProblem

A compiled nominal linear optimization problem.

Returned by: Subclassed by: Used by:

Methods:

  • linprog_args –

    Return coefficients as linprog keyword arguments.

linprog_args 

linprog_args() -> LinprogArgs

Return coefficients as linprog keyword arguments.

Source code in src/vic3_analysis/optimize/base.py 
112
113
114
115
116
117
118
119
120
def linprog_args(self) -> LinprogArgs:
    """Return coefficients as ``linprog`` keyword arguments."""
    return {
        "c": self.objective_vector,
        "A_ub": self.A_ub,
        "b_ub": self.b_ub,
        "A_eq": self.A_eq,
        "b_eq": self.b_eq,
    }

NominalOptimizer 

NominalOptimizer(model: Economy)

Bases: BaseOptimizer[LinearProblem]

Compile and solve nominal-price scenarios as linear programmes.

Methods:

  • compile –

    Compile scenario into an inspectable linear problem.

  • solve –

    Compile and solve scenario, returning its optimal economy state.

  • solve_problem –

    Solve a previously compiled linear problem.

  • import_marginals –

    Return objective sensitivities to increasing imports by good.

Source code in src/vic3_analysis/optimize/base.py 
163
164
165
166
167
168
169
def __init__(self, model: Economy) -> None:
    """Bind the optimizer to *model*."""
    self.model = model
    self.result = None
    self.scenario = None
    self.problem = None
    self._state = None

compile 

compile(scenario: Scenario) -> LinearProblem

Compile scenario into an inspectable linear problem.

Source code in src/vic3_analysis/optimize/nominal.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def compile(self, scenario: Scenario) -> LinearProblem:
    """Compile *scenario* into an inspectable linear problem."""
    if scenario.objective == "gdp_per_capita":
        raise ValueError(
            "NominalOptimizer does not support the nonlinear "
            "gdp_per_capita objective; use MarketOptimizer."
        )
    if scenario.objective == "average_wealth":
        raise ValueError(
            "NominalOptimizer does not support average_wealth; "
            "use EquilibriumOptimizer."
        )
    data = self._compile_common(scenario)
    if scenario.objective == "gdp":
        objective = -data.gdp_vector
    elif scenario.objective == "employment":
        objective = -self.model.employment_vector
    elif scenario.objective == "automation":
        objective = self.model.employment_vector
    elif scenario.objective == "construction_cost":
        objective = self.model.construction_cost_vector
    else:
        raise ValueError(f"Unknown objective: {scenario.objective!r}")
    problem = self._linear_problem(data, objective)
    self._bind_problem(problem)
    return problem

solve 

solve(scenario: Scenario) -> EconomyState

Compile and solve scenario, returning its optimal economy state.

Source code in src/vic3_analysis/optimize/nominal.py 
44
45
46
def solve(self, scenario: Scenario) -> EconomyState:
    """Compile and solve *scenario*, returning its optimal economy state."""
    return self.solve_problem(self.compile(scenario))

solve_problem 

solve_problem(problem: LinearProblem) -> EconomyState

Solve a previously compiled linear problem.

Source code in src/vic3_analysis/optimize/nominal.py 
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def solve_problem(self, problem: LinearProblem) -> EconomyState:
    """Solve a previously compiled linear problem."""
    self._bind_problem(problem)
    result = opt.linprog(
        **problem.linprog_args(),
        bounds=problem.bounds,
        method="highs",
    )
    if not result.success:
        raise ValueError(f"Optimization failed: {result.message}")
    self.result = result
    state = self._state_from_levels(problem, np.asarray(result.x, dtype=np.float64))
    self._state = state
    return state

import_marginals 

import_marginals() -> ndarray | None

Return objective sensitivities to increasing imports by good.

Values align to the goods index and measure the minimized LP objective change per additional imported unit under the goods-balance floor.

Source code in src/vic3_analysis/optimize/nominal.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def import_marginals(self) -> np.ndarray | None:
    """Return objective sensitivities to increasing imports by good.

    Values align to the goods index and measure the minimized LP objective
    change per additional imported unit under the goods-balance floor.
    """
    problem = self.problem
    result = self.result
    if problem is None or result is None:
        return None
    row_slice = problem.inequality_slices.get("min_goods_balance")
    if row_slice is None:
        return None
    ineqlin = getattr(result, "ineqlin", None)
    if ineqlin is None:
        return None
    marginals = getattr(ineqlin, "marginals", None)
    if marginals is None:
        return None
    values = np.asarray(marginals, dtype=np.float64)
    stop = row_slice.stop
    if stop is None or values.shape[0] < stop:
        return None
    return values[row_slice].copy()

MarketProblem dataclass 

MarketProblem(
    economy: Economy,
    scenario: Scenario,
    throughput_multipliers: ndarray,
    input_matrix: ndarray,
    output_matrix: ndarray,
    net_matrix: ndarray,
    gdp_vector: ndarray,
    goods_import: ndarray,
    goods_export: ndarray,
    A_ub: ndarray | None,
    b_ub: ndarray | None,
    A_eq: ndarray | None,
    b_eq: ndarray | None,
    bounds: Bounds,
    inequality_slices: ConstraintSlices,
    equality_slices: ConstraintSlices,
    objective: Callable[[ndarray], float],
    jacobian: Callable[[ndarray], ndarray],
    warm_start: LinearProblem,
)

Bases: _CompiledProblem

A compiled nonlinear market-price optimization problem.

Returned by: Subclassed by: Used by:

MarketOptimizer 

MarketOptimizer(model: Economy)

Bases: BaseOptimizer[MarketProblem]

Compile and solve scenarios with endogenous national market prices.

Methods:

  • compile –

    Compile scenario into an inspectable nonlinear market problem.

  • solve –

    Compile and solve scenario with endogenous market prices.

  • solve_problem –

    Solve a previously compiled market problem.

Source code in src/vic3_analysis/optimize/base.py 
163
164
165
166
167
168
169
def __init__(self, model: Economy) -> None:
    """Bind the optimizer to *model*."""
    self.model = model
    self.result = None
    self.scenario = None
    self.problem = None
    self._state = None

compile 

compile(scenario: Scenario) -> MarketProblem

Compile scenario into an inspectable nonlinear market problem.

Source code in src/vic3_analysis/optimize/market.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def compile(self, scenario: Scenario) -> MarketProblem:
    """Compile *scenario* into an inspectable nonlinear market problem."""
    if scenario.objective not in ("gdp", "gdp_per_capita"):
        raise ValueError(
            "MarketOptimizer requires objective='gdp' or "
            "objective='gdp_per_capita'."
        )
    data = self._compile_common(scenario)
    warm_start_objective = (
        self.model.employment_vector
        if scenario.objective == "gdp_per_capita" and scenario.goods_export
        else -data.gdp_vector
    )
    warm_start = self._extract_fixed_zero_bounds(
        self._linear_problem(data, warm_start_objective)
    )
    employment = self.model.employment_vector
    economy = self.model

    def objective_terms(levels: np.ndarray) -> tuple[float, np.ndarray]:
        value, gradient = self._value_and_gradient(
            economy,
            levels,
            data.input_matrix,
            data.output_matrix,
            data.net_matrix,
            data.goods_import,
            data.goods_export,
            economy.base_prices,
        )
        if scenario.objective == "gdp_per_capita":
            value, gradient = self._per_capita_value_and_gradient(
                value, gradient, levels, employment
            )
            return (
                value * scenario.workforce_ratio,
                gradient * scenario.workforce_ratio,
            )
        return value, gradient

    def objective(levels: np.ndarray) -> float:
        value, _gradient = objective_terms(levels)
        return -value

    def jacobian(levels: np.ndarray) -> np.ndarray:
        _value, gradient = objective_terms(levels)
        return -gradient

    problem = MarketProblem(
        economy=data.economy,
        scenario=scenario,
        throughput_multipliers=data.throughput_multipliers,
        input_matrix=data.input_matrix,
        output_matrix=data.output_matrix,
        net_matrix=data.net_matrix,
        gdp_vector=data.gdp_vector,
        goods_import=data.goods_import,
        goods_export=data.goods_export,
        A_ub=warm_start.A_ub,
        b_ub=warm_start.b_ub,
        A_eq=warm_start.A_eq,
        b_eq=warm_start.b_eq,
        bounds=warm_start.bounds,
        inequality_slices=warm_start.inequality_slices,
        equality_slices=warm_start.equality_slices,
        objective=objective,
        jacobian=jacobian,
        warm_start=warm_start,
    )
    self._bind_problem(problem)
    return problem

solve 

solve(
    scenario: Scenario,
    *,
    x0: ndarray | None = None,
    options: Mapping[str, object] | None = None,
) -> EconomyState

Compile and solve scenario with endogenous market prices.

Source code in src/vic3_analysis/optimize/market.py 
 96
 97
 98
 99
100
101
102
103
104
def solve(
    self,
    scenario: Scenario,
    *,
    x0: np.ndarray | None = None,
    options: Mapping[str, object] | None = None,
) -> EconomyState:
    """Compile and solve *scenario* with endogenous market prices."""
    return self.solve_problem(self.compile(scenario), x0=x0, options=options)

solve_problem 

solve_problem(
    problem: MarketProblem,
    *,
    x0: ndarray | None = None,
    options: Mapping[str, object] | None = None,
) -> EconomyState

Solve a previously compiled market problem.

Source code in src/vic3_analysis/optimize/market.py 
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def solve_problem(
    self,
    problem: MarketProblem,
    *,
    x0: np.ndarray | None = None,
    options: Mapping[str, object] | None = None,
) -> EconomyState:
    """Solve a previously compiled market problem."""
    self._bind_problem(problem)
    warm_start_problem = problem.warm_start
    self._check_feasibility_and_boundedness(warm_start_problem)
    n_buildings = len(self.model.building_index)

    if x0 is None:
        warm_start = opt.linprog(
            **warm_start_problem.linprog_args(),
            bounds=warm_start_problem.bounds,
            method="highs",
        )
        if not warm_start.success:
            raise ValueError(f"Nominal warm-start failed: {warm_start.message}")
        initial = np.asarray(warm_start.x, dtype=np.float64)
        if (
            problem.scenario.objective == "gdp_per_capita"
            and problem.scenario.goods_export
            and float(initial @ self.model.employment_vector) <= 0
        ):
            fallback = replace(
                warm_start_problem, objective_vector=-problem.gdp_vector
            )
            warm_start = opt.linprog(
                **fallback.linprog_args(), bounds=fallback.bounds, method="highs"
            )
            if not warm_start.success:
                raise ValueError(f"Nominal warm-start failed: {warm_start.message}")
            initial = np.asarray(warm_start.x, dtype=np.float64)
    else:
        initial = self._validate_initial_point(
            x0,
            problem.A_ub,
            problem.b_ub,
            problem.A_eq,
            problem.b_eq,
            problem.bounds,
            n_buildings,
        )

    constraints: list[opt.LinearConstraint] = []
    if problem.A_ub is not None and problem.b_ub is not None:
        linear_constraint = cast(Any, opt.LinearConstraint)
        constraints.append(linear_constraint(problem.A_ub, -np.inf, problem.b_ub))
    if problem.A_eq is not None and problem.b_eq is not None:
        linear_constraint = cast(Any, opt.LinearConstraint)
        constraints.append(
            linear_constraint(problem.A_eq, problem.b_eq, problem.b_eq)
        )

    result = opt.minimize(
        problem.objective,
        initial,
        jac=problem.jacobian,
        method="SLSQP",
        bounds=problem.bounds,
        constraints=constraints,
        options=dict(options) if options is not None else None,
    )
    if not result.success:
        raise ValueError(f"Market optimization failed: {result.message}")

    self.result = result
    state = self._state_from_levels(
        problem,
        np.asarray(result.x, dtype=np.float64),
        method="market",
    )
    self._state = state
    return state

EquilibriumProblem dataclass 

EquilibriumProblem(
    economy: Economy,
    scenario: Scenario,
    throughput_multipliers: ndarray,
    input_matrix: ndarray,
    output_matrix: ndarray,
    net_matrix: ndarray,
    gdp_vector: ndarray,
    goods_import: ndarray,
    goods_export: ndarray,
    A_ub: ndarray | None,
    b_ub: ndarray | None,
    A_eq: ndarray | None,
    b_eq: ndarray | None,
    bounds: Bounds,
    inequality_slices: ConstraintSlices,
    equality_slices: ConstraintSlices,
    warm_start: LinearProblem,
    initial_wealth: ndarray,
    max_iterations: int = 1000,
    rtol: float = 1e-09,
    atol: float = 1e-08,
    *,
    subsidize_operating_losses: bool = True,
)

Bases: _CompiledProblem

An immutable equilibrium-search recipe bound to one economy.

Arrays retain full economy index order. initial_wealth is copied and reused for every evaluation, making results independent of earlier trials. warm_start is a nominal GDP LP over the same linear feasible region. Evaluation rebuilds supply-dependent baskets for each building mix.

Returned by: Subclassed by: Used by:

Methods:

  • evaluate –

    Solve household equilibrium at full-index nonnegative building levels.

  • objective –

    Return the negative, unscaled metric after solving equilibrium.

evaluate 

evaluate(levels: ndarray) -> EconomyState

Solve household equilibrium at full-index nonnegative building levels.

Scenario feasibility is checked by the optimizer, separately from this physical evaluation. Invalid model inputs raise ValueError; failed equilibrium iteration raises EquilibriumError with its snapshot.

Source code in src/vic3_analysis/optimize/equilibrium.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def evaluate(self, levels: np.ndarray) -> EconomyState:
    """Solve household equilibrium at full-index nonnegative building levels.

    Scenario feasibility is checked by the optimizer, separately from this
    physical evaluation. Invalid model inputs raise ``ValueError``; failed
    equilibrium iteration raises ``EquilibriumError`` with its snapshot.
    """
    levels = np.asarray(levels, dtype=np.float64)
    if levels.shape != (len(self.bounds),):
        raise ValueError("levels must be aligned to building configurations.")
    if not np.isfinite(levels).all() or (levels < 0).any():
        raise ValueError("levels must contain finite non-negative values.")
    return self.economy.solve(
        levels,
        method="equilibrium",
        goods_import=self.goods_import,
        goods_export=self.goods_export,
        throughput_multipliers=self.throughput_multipliers,
        workforce_ratio=self.scenario.workforce_ratio,
        initial_wealth=self.initial_wealth,
        subsidize_operating_losses=self.subsidize_operating_losses,
        max_iterations=self.max_iterations,
        rtol=self.rtol,
        atol=self.atol,
    )

objective 

objective(levels: ndarray) -> float

Return the negative, unscaled metric after solving equilibrium.

Source code in src/vic3_analysis/optimize/equilibrium.py 
109
110
111
def objective(self, levels: np.ndarray) -> float:
    """Return the negative, unscaled metric after solving equilibrium."""
    return -_metric(self.evaluate(levels), self.scenario.objective)

EquilibriumOptimizer 

EquilibriumOptimizer(model: Economy)

Bases: BaseOptimizer[EquilibriumProblem]

Search continuous building levels using verified household equilibria.

This is a deterministic local heuristic, without a global-optimum or smooth-objective stationarity guarantee. Population varies with employment. Numerical equilibrium permits the existing financial warning flags.

Subclassed by:

Methods:

  • summary –

    Extend the solution report with search and equilibrium diagnostics.

  • compile –

    Compile shared constraints and freeze the equilibrium evaluation seed.

  • solve –

    Compile and search a scenario, returning the best verified equilibrium.

  • solve_problem –

    Search a compiled problem over a bounded, feasible linear region.

Source code in src/vic3_analysis/optimize/base.py 
163
164
165
166
167
168
169
def __init__(self, model: Economy) -> None:
    """Bind the optimizer to *model*."""
    self.model = model
    self.result = None
    self.scenario = None
    self.problem = None
    self._state = None

summary 

summary() -> DataFrame

Extend the solution report with search and equilibrium diagnostics.

Evaluation counts exclude final verification, matching result. Financial warning counts do not imply equilibrium non-convergence.

Source code in src/vic3_analysis/optimize/equilibrium.py 
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
def summary(self) -> pd.DataFrame:
    """Extend the solution report with search and equilibrium diagnostics.

    Evaluation counts exclude final verification, matching ``result``.
    Financial warning counts do not imply equilibrium non-convergence.
    """
    report = super().summary()
    diagnostics = self._require_state().equilibrium_diagnostics
    result = self.result
    if diagnostics is None or result is None:
        raise ValueError("Optimizer has no returned equilibrium to summarize.")
    return report.assign(
        equilibrium_evaluations=int(result.equilibrium_evaluations),
        failed_equilibrium_evaluations=sum(result.equilibrium_failures.values()),
        equilibrium_status=diagnostics.status,
        equilibrium_iterations=diagnostics.iterations,
        loss_buildings_count=len(diagnostics.loss_buildings),
        unaffordable_floor_cells_count=len(diagnostics.unaffordable_floor_cells),
    )

compile 

compile(
    scenario: Scenario,
    *,
    initial_wealth: ndarray | None = None,
    subsidize_operating_losses: bool = True,
    max_iterations: int = 1000,
    rtol: float = 1e-09,
    atol: float = 1e-08,
) -> EquilibriumProblem

Compile shared constraints and freeze the equilibrium evaluation seed.

Supports gdp, gdp_per_capita, and average_wealth. Household consumption is endogenous; imports and exports are fixed external orders. Goods-balance floors include trade and exclude household consumption. subsidize_operating_losses funds losses through proportional taxes on cooperative and dependent income, without guaranteeing payroll.

Source code in src/vic3_analysis/optimize/equilibrium.py 
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
def compile(
    self,
    scenario: Scenario,
    *,
    initial_wealth: np.ndarray | None = None,
    subsidize_operating_losses: bool = True,
    max_iterations: int = 1000,
    rtol: float = 1e-9,
    atol: float = 1e-8,
) -> EquilibriumProblem:
    """Compile shared constraints and freeze the equilibrium evaluation seed.

    Supports ``gdp``, ``gdp_per_capita``, and ``average_wealth``. Household
    consumption is endogenous; imports and exports are fixed external orders.
    Goods-balance floors include trade and exclude household consumption.
    ``subsidize_operating_losses`` funds losses through proportional taxes
    on cooperative and dependent income, without guaranteeing payroll.
    """
    if scenario.objective not in ("gdp", "gdp_per_capita", "average_wealth"):
        raise ValueError(
            "EquilibriumOptimizer requires objective='gdp', "
            "'gdp_per_capita', or 'average_wealth'."
        )
    data = self._compile_common(scenario)
    warm_start = self._extract_fixed_zero_bounds(
        self._linear_problem(data, -data.gdp_vector)
    )
    # Exact zero caps also imply fixed variables. Keep the named inequality
    # rows for inspection; only the search coordinates are reduced.
    if warm_start.A_ub is not None and warm_start.b_ub is not None:
        bounds = list(warm_start.bounds)
        for row, rhs in zip(warm_start.A_ub, warm_start.b_ub, strict=True):
            if rhs == 0 and (row >= 0).all():
                for index in np.flatnonzero(row > 0):
                    bounds[int(index)] = (0.0, 0.0)
        warm_start = replace(warm_start, bounds=tuple(bounds))
    if initial_wealth is None:
        initial_wealth = np.broadcast_to(
            self.model.pop_wealth_init,
            (len(self.model.building_index), len(self.model.pop_index)),
        )
    problem = EquilibriumProblem(
        economy=data.economy,
        scenario=scenario,
        throughput_multipliers=data.throughput_multipliers,
        input_matrix=data.input_matrix,
        output_matrix=data.output_matrix,
        net_matrix=data.net_matrix,
        gdp_vector=data.gdp_vector,
        goods_import=data.goods_import,
        goods_export=data.goods_export,
        A_ub=warm_start.A_ub,
        b_ub=warm_start.b_ub,
        A_eq=warm_start.A_eq,
        b_eq=warm_start.b_eq,
        bounds=warm_start.bounds,
        inequality_slices=warm_start.inequality_slices,
        equality_slices=warm_start.equality_slices,
        warm_start=warm_start,
        initial_wealth=initial_wealth,
        subsidize_operating_losses=subsidize_operating_losses,
        max_iterations=max_iterations,
        rtol=rtol,
        atol=atol,
    )
    self._bind_problem(problem)
    return problem

solve 

solve(
    scenario: Scenario,
    *,
    x0: ndarray | None = None,
    options: Mapping[str, object] | None = None,
    log_every: int = 100,
    initial_wealth: ndarray | None = None,
    subsidize_operating_losses: bool = True,
    max_iterations: int = 1000,
    rtol: float = 1e-09,
    atol: float = 1e-08,
) -> EconomyState

Compile and search a scenario, returning the best verified equilibrium.

options accepts COBYLA's rhobeg, tol, catol, disp, and f_target plus a hard maxiter trial limit (default 20,000). max_iterations separately limits each household equilibrium solve. subsidize_operating_losses enables income-tax-funded operating losses; trials needing more than all household income are rejected. log_every sets the trial interval for INFO progress summaries when this module's logger is enabled; it must be a positive integer. See solve_problem for search termination and failure diagnostics.

Source code in src/vic3_analysis/optimize/equilibrium.py 
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
def solve(
    self,
    scenario: Scenario,
    *,
    x0: np.ndarray | None = None,
    options: Mapping[str, object] | None = None,
    log_every: int = 100,
    initial_wealth: np.ndarray | None = None,
    subsidize_operating_losses: bool = True,
    max_iterations: int = 1000,
    rtol: float = 1e-9,
    atol: float = 1e-8,
) -> EconomyState:
    """Compile and search a scenario, returning the best verified equilibrium.

    ``options`` accepts COBYLA's ``rhobeg``, ``tol``, ``catol``, ``disp``,
    and ``f_target`` plus a hard ``maxiter`` trial limit (default 20,000).
    ``max_iterations`` separately limits each household equilibrium solve.
    ``subsidize_operating_losses`` enables income-tax-funded operating losses;
    trials needing more than all household income are rejected.
    ``log_every`` sets the trial interval for INFO progress summaries when
    this module's logger is enabled; it must be a positive integer.
    See ``solve_problem`` for search termination and failure diagnostics.
    """
    problem = self.compile(
        scenario,
        initial_wealth=initial_wealth,
        subsidize_operating_losses=subsidize_operating_losses,
        max_iterations=max_iterations,
        rtol=rtol,
        atol=atol,
    )
    return self.solve_problem(problem, x0=x0, options=options, log_every=log_every)

solve_problem 

solve_problem(
    problem: EquilibriumProblem,
    *,
    x0: ndarray | None = None,
    options: Mapping[str, object] | None = None,
    log_every: int = 100,
) -> EconomyState

Search a compiled problem over a bounded, feasible linear region.

Return the best feasible equilibrium even on budget exhaustion; inspect result.success and result.message for search termination. result.x and unscaled result.fun describe the returned state. result.solver_result retains raw COBYLA output (or None if it was skipped or interrupted by the hard budget). nfev counts distinct consecutive trials, including rejected points; equilibrium_evaluations counts inner solve attempts. The initial trial counts against the budget; final verification does not. equilibrium_failures counts failures by status and last_equilibrium_error retains the last failed inner solve.

Enable vic3_analysis.optimize.equilibrium at INFO for lifecycle events and summaries every log_every trials (default 100, positive integer). DEBUG additionally reports uncached trials and incumbents. Logging uses actual maximizing metrics in explicit units, independently of native COBYLA disp. Library code never configures log handlers.

Raises:

  • ValueError –

    For invalid inputs, infeasible/unbounded regions, or no feasible converged equilibrium (with diagnostics on result).

  • EquilibriumError –

    If fresh verification of the selected state fails.

Source code in src/vic3_analysis/optimize/equilibrium.py 
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
def solve_problem(
    self,
    problem: EquilibriumProblem,
    *,
    x0: np.ndarray | None = None,
    options: Mapping[str, object] | None = None,
    log_every: int = 100,
) -> EconomyState:
    """Search a compiled problem over a bounded, feasible linear region.

    Return the best feasible equilibrium even on budget exhaustion; inspect
    ``result.success`` and ``result.message`` for search termination.
    ``result.x`` and unscaled ``result.fun`` describe the returned state.
    ``result.solver_result`` retains raw COBYLA output (or ``None`` if it
    was skipped or interrupted by the hard budget). ``nfev`` counts
    distinct consecutive trials, including
    rejected points; ``equilibrium_evaluations`` counts inner solve attempts.
    The initial trial counts against the budget; final verification does not.
    ``equilibrium_failures`` counts failures by status and
    ``last_equilibrium_error`` retains the last failed inner solve.

    Enable ``vic3_analysis.optimize.equilibrium`` at INFO for lifecycle
    events and summaries every ``log_every`` trials (default 100, positive
    integer). DEBUG additionally reports uncached trials and incumbents.
    Logging uses actual maximizing metrics in explicit units, independently
    of native COBYLA ``disp``. Library code never configures log handlers.

    Raises:
        ValueError: For invalid inputs, infeasible/unbounded regions, or no
            feasible converged equilibrium (with diagnostics on ``result``).
        EquilibriumError: If fresh verification of the selected state fails.
    """
    self._bind_problem(problem)
    solver_options, budget = _options(options)
    if (
        isinstance(log_every, (bool, np.bool_))
        or not isinstance(log_every, (int, np.integer))
        or log_every <= 0
    ):
        raise ValueError("log_every must be a positive integer.")
    search = _Search(problem, budget, log_every=int(log_every), phase="preflight")
    search.report("started")
    try:
        return self._solve_problem(problem, x0, solver_options, search)
    except Exception as error:
        search.report("finished", status="failed", message=str(error))
        raise

SupplyChainEquilibriumProblem dataclass 

SupplyChainEquilibriumProblem(
    economy: Economy,
    scenario: Scenario,
    throughput_multipliers: ndarray,
    input_matrix: ndarray,
    output_matrix: ndarray,
    net_matrix: ndarray,
    gdp_vector: ndarray,
    goods_import: ndarray,
    goods_export: ndarray,
    A_ub: ndarray | None,
    b_ub: ndarray | None,
    A_eq: ndarray | None,
    b_eq: ndarray | None,
    bounds: Bounds,
    inequality_slices: ConstraintSlices,
    equality_slices: ConstraintSlices,
    warm_start: LinearProblem,
    initial_wealth: ndarray,
    max_iterations: int = 1000,
    rtol: float = 1e-09,
    atol: float = 1e-08,
    *,
    subsidize_operating_losses: bool = True,
    chain_index: tuple[str, ...],
    chain_matrix: ndarray,
    building_problem: EquilibriumProblem,
    _owner: object = object(),
)

Bases: EquilibriumProblem

Chain-coordinate constraints with building-coordinate economic accounting.

bounds, constraint columns, and warm_start decision coefficients follow chain_index. Production matrices and wealth retain the original economy indexes. Physical evaluation delegates to building_problem.

Returned by:

Methods:

  • to_building_levels –

    Map chain coordinates using this problem's immutable coefficients.

  • evaluate –

    Solve equilibrium after transforming a nonnegative chain vector.

to_building_levels 

to_building_levels(chain_levels: ndarray) -> ndarray

Map chain coordinates using this problem's immutable coefficients.

Source code in src/vic3_analysis/optimize/supply_chain.py 
79
80
81
def to_building_levels(self, chain_levels: np.ndarray) -> np.ndarray:
    """Map chain coordinates using this problem's immutable coefficients."""
    return _to_building_levels(self.chain_matrix, chain_levels)

evaluate 

evaluate(levels: ndarray) -> EconomyState

Solve equilibrium after transforming a nonnegative chain vector.

Source code in src/vic3_analysis/optimize/supply_chain.py 
83
84
85
def evaluate(self, levels: np.ndarray) -> EconomyState:
    """Solve equilibrium after transforming a nonnegative chain vector."""
    return self.building_problem.evaluate(self.to_building_levels(levels))

SupplyChainEquilibriumOptimizer 

SupplyChainEquilibriumOptimizer(
    economy: Economy, df_building_levels: DataFrame
)

Bases: EquilibriumOptimizer

Optimize nonnegative chain levels using the existing equilibrium search.

x0 and result.x contain chain levels. Returned states and result.building_levels contain the recipe matrix times result.x. Fixed recipes restrict the search space; they do not guarantee financial viability.

DataFrame rows are aligned by configuration key to the economy; all-NaN columns are skipped with a warning. Shared configurations receive the sum of chain contributions, without normalization of coefficients.

Methods:

  • summary –

    Add the recipe count to the report of the verified building state.

  • validate_building_levels –

    Validate, align, and copy recipes, warning when failed columns are skipped.

  • to_building_levels –

    Map a finite nonnegative chain vector to building configurations.

  • compile –

    Compile building constraints and the nominal warm start into chains.

  • solve_problem –

    Search chain coordinates and verify the original building constraints.

Attributes:

  • chain_index (tuple[str, ...]) –

    Return surviving recipe labels in chain-vector order.

  • df_building_levels (DataFrame) –

    Return an independent copy of the aligned production-unit-by-chain matrix.

Source code in src/vic3_analysis/optimize/supply_chain.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def __init__(self, economy: Economy, df_building_levels: pd.DataFrame) -> None:
    """Snapshot a production-unit-by-chain DataFrame for an economy.

    DataFrame rows are aligned by configuration key to the economy;
    all-NaN columns are skipped with a warning. Shared configurations receive
    the sum of chain contributions, without normalization of coefficients.
    """
    super().__init__(economy)
    frame = self.validate_building_levels(
        df_building_levels, building_index=tuple(economy.building_index)
    )
    self._chain_index = tuple(frame.columns)
    self._chain_matrix = _readonly_array(frame.to_numpy(dtype=np.float64))
    self._owner = object()

chain_index property 

chain_index: tuple[str, ...]

Return surviving recipe labels in chain-vector order.

df_building_levels property 

df_building_levels: DataFrame

Return an independent copy of the aligned production-unit-by-chain matrix.

summary 

summary() -> DataFrame

Add the recipe count to the report of the verified building state.

chain_count counts recipe columns; configuration_count counts building configurations. Economy metrics use mapped building levels.

Source code in src/vic3_analysis/optimize/supply_chain.py 
115
116
117
118
119
120
121
def summary(self) -> pd.DataFrame:
    """Add the recipe count to the report of the verified building state.

    ``chain_count`` counts recipe columns; ``configuration_count`` counts
    building configurations. Economy metrics use mapped building levels.
    """
    return super().summary().assign(chain_count=len(self.chain_index))

validate_building_levels staticmethod 

validate_building_levels(
    frame: DataFrame, *, building_index: tuple[str, ...]
) -> DataFrame

Validate, align, and copy recipes, warning when failed columns are skipped.

Rows must match building_index exactly; column order is preserved. All-NaN columns are omitted. Each remaining recipe must have finite, nonnegative coefficients and at least one positive coefficient. Returned float coefficients retain the caller's units without normalization.

Source code in src/vic3_analysis/optimize/supply_chain.py 
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
@staticmethod
def validate_building_levels(
    frame: pd.DataFrame,
    *,
    building_index: tuple[str, ...],
) -> pd.DataFrame:
    """Validate, align, and copy recipes, warning when failed columns are skipped.

    Rows must match ``building_index`` exactly; column order is preserved.
    All-NaN columns are omitted. Each remaining recipe must have finite,
    nonnegative coefficients and at least one positive coefficient. Returned
    float coefficients retain the caller's units without normalization.
    """
    buildings = tuple(building_index)
    rows = tuple(frame.index)
    columns = tuple(frame.columns)
    for name, labels in (
        ("building_index", buildings),
        ("DataFrame index", rows),
        ("DataFrame columns", columns),
    ):
        if any(not isinstance(key, str) or not key for key in labels):
            raise ValueError(f"{name} must contain nonempty string labels.")
        if len(set(labels)) != len(labels):
            raise ValueError(f"{name} must contain unique labels.")
    if not buildings or set(rows) != set(buildings):
        raise ValueError("DataFrame index must match building_index exactly.")
    failed = frame.isna().all(axis=0)
    skipped = [key for key, missing in zip(columns, failed, strict=True) if missing]
    if skipped:
        warnings.warn(
            "Skipping failed supply-chain recipes containing only NaN: "
            + ", ".join(skipped),
            UserWarning,
            stacklevel=2,
        )
    chains = tuple(key for key in columns if key not in skipped)
    if not chains:
        raise ValueError("No usable supply-chain recipes remain.")
    try:
        matrix = frame.loc[list(buildings), list(chains)].to_numpy(
            dtype=np.float64, copy=True
        )
    except (TypeError, ValueError) as exc:
        raise ValueError("matrix must contain numeric coefficients.") from exc
    if not np.isfinite(matrix).all() or (matrix < 0).any():
        raise ValueError("matrix must contain finite non-negative coefficients.")
    if not (matrix > 0).any(axis=0).all():
        raise ValueError("Every chain must have at least one nonzero coefficient.")
    return pd.DataFrame(
        matrix, index=pd.Index(buildings, name="config"), columns=chains
    )

to_building_levels 

to_building_levels(chain_levels: ndarray) -> ndarray

Map a finite nonnegative chain vector to building configurations.

Source code in src/vic3_analysis/optimize/supply_chain.py 
190
191
192
def to_building_levels(self, chain_levels: np.ndarray) -> np.ndarray:
    """Map a finite nonnegative chain vector to building configurations."""
    return _to_building_levels(self._chain_matrix, chain_levels)

compile 

compile(
    scenario: Scenario,
    *,
    initial_wealth: ndarray | None = None,
    subsidize_operating_losses: bool = True,
    max_iterations: int = 1000,
    rtol: float = 1e-09,
    atol: float = 1e-08,
) -> SupplyChainEquilibriumProblem

Compile building constraints and the nominal warm start into chains.

Source code in src/vic3_analysis/optimize/supply_chain.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
def compile(
    self,
    scenario: Scenario,
    *,
    initial_wealth: np.ndarray | None = None,
    subsidize_operating_losses: bool = True,
    max_iterations: int = 1000,
    rtol: float = 1e-9,
    atol: float = 1e-8,
) -> SupplyChainEquilibriumProblem:
    """Compile building constraints and the nominal warm start into chains."""
    building = super().compile(
        scenario,
        initial_wealth=initial_wealth,
        subsidize_operating_losses=subsidize_operating_losses,
        max_iterations=max_iterations,
        rtol=rtol,
        atol=atol,
    )
    matrix = self._chain_matrix
    A_ub = None if building.A_ub is None else building.A_ub @ matrix
    A_eq = None if building.A_eq is None else building.A_eq @ matrix
    b_ub = building.b_ub
    slices = dict(building.inequality_slices)
    bounds: Bounds = tuple((0.0, None) for _ in self.chain_index)
    blocked = np.zeros(len(bounds), dtype=bool)
    rows: list[np.ndarray] = []
    rhs: list[float] = []
    for row, (lower, upper) in zip(matrix, building.bounds, strict=True):
        if upper == 0:
            blocked |= row > 0
        elif upper is not None:
            rows.append(row)
            rhs.append(upper)
        if lower > 0:
            rows.append(-row)
            rhs.append(-lower)
    if rows:
        start = 0 if b_ub is None else len(b_ub)
        extra = np.stack(rows)
        A_ub = extra if A_ub is None else np.vstack((A_ub, extra))
        b_ub = np.asarray(rhs) if b_ub is None else np.concatenate((b_ub, rhs))
        slices["building_bounds"] = slice(start, start + len(rhs))
    bounds = tuple((0.0, 0.0 if used else None) for used in blocked)
    warm = self._extract_fixed_zero_bounds(
        replace(
            building.warm_start,
            A_ub=A_ub,
            b_ub=b_ub,
            A_eq=A_eq,
            bounds=bounds,
            inequality_slices=slices,
            objective_vector=matrix.T @ building.warm_start.objective_vector,
        )
    )
    problem = SupplyChainEquilibriumProblem(
        economy=building.economy,
        scenario=building.scenario,
        throughput_multipliers=building.throughput_multipliers,
        input_matrix=building.input_matrix,
        output_matrix=building.output_matrix,
        net_matrix=building.net_matrix,
        gdp_vector=building.gdp_vector,
        goods_import=building.goods_import,
        goods_export=building.goods_export,
        A_ub=warm.A_ub,
        b_ub=warm.b_ub,
        A_eq=warm.A_eq,
        b_eq=warm.b_eq,
        bounds=warm.bounds,
        inequality_slices=warm.inequality_slices,
        equality_slices=warm.equality_slices,
        warm_start=warm,
        initial_wealth=building.initial_wealth,
        subsidize_operating_losses=building.subsidize_operating_losses,
        max_iterations=building.max_iterations,
        rtol=building.rtol,
        atol=building.atol,
        chain_index=self.chain_index,
        chain_matrix=matrix,
        _owner=self._owner,
        building_problem=building,
    )
    self._bind_problem(problem)
    return problem

solve_problem 

solve_problem(
    problem: EquilibriumProblem,
    *,
    x0: ndarray | None = None,
    options: Mapping[str, object] | None = None,
    log_every: int = 100,
) -> EconomyState

Search chain coordinates and verify the original building constraints.

Source code in src/vic3_analysis/optimize/supply_chain.py 
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
def solve_problem(
    self,
    problem: EquilibriumProblem,
    *,
    x0: np.ndarray | None = None,
    options: Mapping[str, object] | None = None,
    log_every: int = 100,
) -> EconomyState:
    """Search chain coordinates and verify the original building constraints."""
    if not isinstance(problem, SupplyChainEquilibriumProblem):
        raise TypeError("Expected a SupplyChainEquilibriumProblem.")
    if (
        problem._owner is not self._owner
        or problem.chain_index != self.chain_index
        or not np.array_equal(problem.chain_matrix, self._chain_matrix)
    ):
        raise ValueError(
            "Compiled problem belongs to a different supply-chain optimizer or recipes."
        )
    try:
        state = super().solve_problem(
            problem, x0=x0, options=options, log_every=log_every
        )
    finally:
        if self.result is not None:
            self.result.building_levels = (
                None
                if self.result.x is None
                else problem.to_building_levels(self.result.x)
            )
    if not _is_feasible(problem.building_problem, state.building_levels):
        self._state = None
        if self.result is not None:
            self.result.success = False
            self.result.message = (
                "Supply-chain final building feasibility check failed."
            )
        raise ValueError("Supply-chain final building feasibility check failed.")
    return state