"""
Profit maximization oracles with Cobb-Douglas production functions.
This module implements oracles for profit maximization where the production
function is of the Cobb-Douglas type: q = A ⋅ y₁^α ⋅ y₂^β.
Three oracle variants:
- ProfitOracle: Standard continuous optimization over log-space variables.
- ProfitRbOracle: Robust variant handling parameter uncertainties in
elasticities, prices, and production limits.
- ProfitQOracle: Discrete (integer-quantized) variant that rounds
continuous solutions and manages infeasibility via retry mechanism.
All oracles return cutting planes (subgradient + violation) that the
cutting-plane algorithm uses to iteratively narrow the search space.
"""
import copy
import math
from typing import Optional, Tuple
import numpy as np
from ellalgo.cutting_plane import OracleOptim, OracleOptimQ
Arr = np.ndarray
Cut = Tuple[Arr, float]
[docs]
class ProfitOracle(OracleOptim):
"""
Oracle for a profit maximization problem with a Cobb-Douglas production function.
This class implements the `OracleOptim` interface for a specific profit
maximization problem. The production function is of the Cobb-Douglas type,
which is widely used in economics to represent the relationship between
production inputs and the amount of output.
The optimization problem is to maximize the profit, which is the difference
between the revenue from selling the product and the cost of the inputs.
The problem is subject to a constraint on one of the inputs.
The `assess_optim` method is the core of the oracle. It takes a candidate
solution (a vector of input quantities) and the current best profit, and
it returns a cutting plane that helps to narrow down the search for the
optimal solution.
"""
idx: int = -1 # Index for round-robin constraint checking
log_Cobb: float # Log value of Cobb-Douglas production
q: Arr # Intermediate calculation of price_out * exp(y)
vy: float # Total variable cost v₁y₁ + v₂y₂
log_pA: float # log(p*A) precomputed value
log_k: float # log(k) constraint value
price_out: Arr # Output prices [v₁, v₂]
elasticities: Arr # Elasticity parameters [α, β]
def __init__(
self, params: Tuple[float, float, float], elasticities: Arr, price_out: Arr
) -> None:
"""Initialize profit maximization oracle with problem parameters.
Parameters:
:param params: Tuple containing:
- unit_price (p): Price per output unit
- scale (A): Production scale factor
- limit (k): Upper bound for x₁
:param elasticities: Array [α, β] of output elasticities
:param price_out: Array [v₁, v₂] of input prices
Mathematical precomputations:
- log_pA = log(p*A) simplifies subsequent exponential calculations
- log_k = log(k) enables log-space constraint checking
"""
unit_price, scale, limit = params
self.log_pA = math.log(unit_price * scale)
self.log_k = math.log(limit)
self.price_out = price_out
self.elasticities = elasticities
self.fns = (self.fn1, self.fn2) # Constraint functions
self.grads = (self.grad1, self.grad2) # Gradient functions
[docs]
def fn1(self, x: Arr, _: float) -> float:
"""Constraint function for y₁ ≤ k (in log-space).
Args:
x: Log-scale input vector [log(y₁), log(y₂)]
Returns:
Constraint violation measure: x[0] - log(k)
Positive values indicate constraint violation
"""
return x[0] - self.log_k # log(y₁) ≤ log(k) → y₁ ≤ k
[docs]
def fn2(self, x: Arr, gamma: float) -> float:
"""Optimality condition function for profit maximization.
Computes:
- Cobb-Douglas value in log-space: log(pA) + αlog(y₁) + βlog(y₂)
- Variable costs: v₁y₁ + v₂y₂
- Optimality gap: log(γ + vy) - log_Cobb
Args:
x: Log-scale input vector
gamma: Current best profit estimate
Updates intermediate values used in gradient calculations
"""
self.log_Cobb = self.log_pA + self.elasticities.dot(x)
self.q = self.price_out * np.exp(x) # v₁y₁, v₂y₂
self.vy = self.q[0] + self.q[1] # Total variable cost
return math.log(gamma + self.vy) - self.log_Cobb
[docs]
def grad1(self, _: float) -> Arr:
"""Gradient for y₁ ≤ k constraint.
Returns:
Gradient vector [1, 0] since ∂(x₀ - log_k)/∂x = (1, 0)
"""
return np.array([1.0, 0.0])
[docs]
def grad2(self, gamma: float) -> Arr:
"""Gradient of optimality condition function.
Computes:
∇f = [v₁y₁/(γ+vy) - α, v₂y₂/(γ+vy) - β]
Args:
gamma: Current profit estimate used in denominator
Uses precomputed q (v₁y₁, v₂y₂) from last fn2 call
"""
return self.q / (gamma + self.vy) - self.elasticities
[docs]
def assess_feas(self, xc: Arr, gamma: float) -> Optional[Cut]:
"""Feasibility assessment using round-robin constraint checking.
Implements:
- Alternates between checking y₁ constraint (fn1) and optimality (fn2)
- Returns first violated constraint found
Args:
xc: Current solution point in log-space
gamma: Current best profit estimate
Returns:
Cut (gradient, violation) if constraint violated
None if all constraints satisfied
"""
for _ in [0, 1]:
self.idx += 1
if self.idx == 2:
self.idx = 0 # Round-robin reset
if (fj := self.fns[self.idx](xc, gamma)) > 0:
return self.grads[self.idx](gamma), fj
return None
[docs]
def assess_optim(self, xc: Arr, gamma: float) -> Tuple[Cut, Optional[float]]:
"""
Assess the optimality of a candidate solution `xc`.
This method is the core of the `ProfitOracle`. It takes a candidate
solution `xc` and the current best profit `gamma`, and it returns a
cutting plane that helps to narrow down the search for the optimal
solution.
The method first checks if the solution is feasible. If not, it returns
a feasibility cut. If the solution is feasible, it calculates the
profit at `xc` and generates an optimality cut.
Args:
xc (Arr): The candidate solution vector (in log-space).
gamma (float): The current best profit.
Returns:
Tuple[Cut, Optional[float]]: A tuple containing the cutting plane
`(g, beta)` and the updated profit `gamma_new`. If the solution is
infeasible, `gamma_new` is `None`.
"""
cut = self.assess_feas(xc, gamma)
if cut is not None:
return cut, None
# Calculate new profit estimate: pA x^α - vy
gamma = np.exp(self.log_Cobb) - self.vy
grad = self.q / (gamma + self.vy) - self.elasticities
return (grad, 0.0), gamma
[docs]
class ProfitRbOracle(OracleOptim):
"""Robust profit oracle handling parameter uncertainty.
Implements robust optimization version from [Aliabadi and Salahi, 2013]
considering uncertainties in:
- Elasticity parameters (α, β)
- Price parameters (p, v)
- Production limit (k)
Uses interval-based uncertainty sets for robust constraint satisfaction.
"""
def __init__(
self,
params: Tuple[float, float, float],
elasticities: Arr,
price_out: Arr,
vparams: Tuple[float, float, float, float, float],
) -> None:
"""Initialize robust oracle with uncertainty parameters.
Parameters:
:param vparams: Uncertainty parameters tuple (ε₁, ε₂, ε₃, ε₄, ε₅) representing:
- ε₁, ε₂: Elasticity uncertainties
- ε₃: Price uncertainty
- ε₄: Production limit uncertainty
- ε₅: Input price uncertainty
Constructs worst-case scenario parameters for robust optimization.
"""
e1, e2, e3, e4, e5 = vparams
self.elasticities = elasticities
self.uie = [e1, e2] # Elasticity uncertainties
unit_price, scale, limit = params
# Construct robust parameters:
params_rb = (
unit_price - e3, # Worst-case price decrease
scale,
limit - e4, # Worst-case capacity reduction
)
self.omega = ProfitOracle(
params_rb,
elasticities,
price_out + np.array([e5, e5]), # Worst-case input price increase
)
[docs]
def assess_optim(self, xc: Arr, gamma: float) -> Tuple[Cut, Optional[float]]:
"""Robust optimization assessment accounting for parameter uncertainties.
Adjusts elasticities based on direction of uncertainty impact:
- Decreases effective α, β when y > 0 (conservative adjustment)
- Increases effective α, β when y ≤ 0 (worst-case scenario)
"""
a_rb = copy.copy(self.elasticities)
for i in [0, 1]:
a_rb[i] += -self.uie[i] if xc[i] > 0.0 else self.uie[i]
self.omega.elasticities = a_rb
return self.omega.assess_optim(xc, gamma)
[docs]
class ProfitQOracle(OracleOptimQ):
"""Discrete profit oracle for integer input quantities.
Solves mixed-integer version of the profit maximization problem:
max p(A y₁^α y₂^β) − v₁y₁ − v₂y₂
s.t. x₁ ≤ k, x ∈ ℕ²
Uses continuous relaxation followed by rounding to nearest integer,
with mechanisms to handle infeasible integer solutions.
"""
xd: np.ndarray # Discrete candidate solution in log-space
def __init__(
self, params: Tuple[float, float, float], elasticities: Arr, price_out: Arr
) -> None:
"""Initialize discrete oracle with underlying continuous oracle."""
self.omega = ProfitOracle(params, elasticities, price_out)
self.xd = np.array([0.0, 0.0]) # Initial discrete solution
[docs]
def assess_optim_q(
self, xc: np.ndarray, gamma: float, retry: bool
) -> Tuple[Tuple[np.ndarray, float], np.ndarray, Optional[float], bool]:
"""Discrete optimization assessment with rounding mechanism.
Workflow:
1. First try continuous solution (retry=False)
2. If infeasible, return feasibility cut
3. If feasible, round to nearest integer and check optimality
4. On retry (retry=True), check rounded solution optimality
Returns:
Tuple containing:
- Cut information
- Evaluation point (continuous or rounded)
- Updated gamma value
- Retry flag for integer feasibility check
"""
if not retry:
# First attempt with continuous solution
if cut := self.omega.assess_feas(xc, gamma):
return cut, xc, None, True
# Round to nearest integer (with 0 → 1 protection)
yd = np.round(np.exp(xc))
yd[yd == 0] = 1.0
self.xd = np.log(yd)
# Check optimality of discrete solution
(grad, beta), gamma_new = self.omega.assess_optim(self.xd, gamma)
beta += grad.dot(self.xd - xc) # Adjust for rounding difference
return (grad, beta), self.xd, gamma_new, not retry