"""
Profit Oracle
This code defines several classes that implement oracles for profit
maximization problems. An oracle, in this context, is a function that helps
solve optimization problems by providing information about the feasibility and
optimality of potential solutions.
The main class, ProfitOracle, is designed to solve a specific type of profit
maximization problem. It takes as input parameters related to production (like
price, scale, and limits) and output elasticities. The goal is to find the
optimal input quantities that maximize profit, given certain constraints.
The ProfitOracle class has methods to assess the feasibility of a solution
(assess_feas) and to find the optimal solution (assess_optim). These methods
take as input a vector y (representing input quantities in log scale) and a
gamma value (representing the current best solution). They output "cuts", which
are linear constraints that help narrow down the search for the optimal
solution.
The code also includes two variations of the profit oracle:
1. ProfitRbOracle: This is a robust version of the profit oracle that can handle some uncertainty in the input parameters.
2. ProfitQOracle: This version deals with discrete (integer) input quantities, as opposed to continuous ones.
The main logic flow in these classes involves calculating various economic
functions (like Cobb-Douglas production functions) and their gradients. The code
uses these calculations to determine if a given solution is feasible and to
guide the search towards the optimal solution.
The output of these oracles is typically a "cut" (a linear constraint) and
possibly an updated best solution (gamma). These outputs are used by an
external optimization algorithm (not shown in this code) to iteratively improve
the solution until the optimal one is found.
For beginners, it's important to understand that this code is implementing
mathematical optimization techniques. While the details might be complex, the
basic idea is to efficiently search for the best solution to a profit
maximization problem, given certain constraints and economic relationships.
"""
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