pyomo

What this skill does

# Pyomo - Mathematical Optimization Modeling

Pyomo allows you to define optimization problems using a natural mathematical syntax (Sets, Parameters, Variables, Constraints). It decouples the model from the solver, allowing the same model to be solved by different engines without code changes.

## FIRST: Verify Prerequisites

```bash
pip install pyomo
# Also install a solver (e.g., GLPK for linear/integer problems)
# Conda: conda install -c conda-forge glpk ipopt
```

## When to Use

- **Strategic Planning**: Long-term resource allocation or investment planning.
- **Process Engineering**: Optimizing chemical plants or refinery operations (Non-linear).
- **Energy Systems**: Power grid dispatch and unit commitment problems.
- **Supply Chain Optimization**: Multi-period, multi-commodity flow problems.
- **Non-Linear Programming (NLP)**: When your constraints or objectives involve smooth curves (e.g., x², log(x)).
- **Stochastic Programming**: Modeling uncertainty in optimization.
- **Custom Solver Integration**: When you need to use specific solvers like IPOPT, SCIP, or Baron.

## Reference Documentation

**Official docs**: http://www.pyomo.org/  
**GitHub**: https://github.com/Pyomo/pyomo  
**Search patterns**: `pyo.ConcreteModel`, `pyo.Constraint`, `pyo.Objective`, `pyo.SolverFactory`

## Core Principles

### Concrete vs. Abstract Models

- **ConcreteModel**: Data is defined at the time the model is built (most common in Python/Data Science).
- **AbstractModel**: The structure is defined first, and data is loaded later (standard for large-scale industrial models).

### Components

- **Var**: Unknowns the solver needs to find.
- **Set/Param**: Data that defines the problem instance.
- **Objective**: The function to minimize or maximize.
- **Constraint**: Rules the variables must follow.

### Solvers

Pyomo does not have its own solver. It requires external solvers (like glpk for LP/MIP or ipopt for NLP) installed on the system.

## Quick Reference

### Installation

```bash
pip install pyomo
# Also install a solver (e.g., GLPK for linear/integer problems)
# Conda: conda install -c conda-forge glpk ipopt
```

### Standard Imports

```python
import pyomo.environ as pyo
from pyomo.opt import SolverFactory
```

### Basic Pattern - Simple Linear Model (Concrete)

```python
import pyomo.environ as pyo

# 1. Create Model
model = pyo.ConcreteModel()

# 2. Define Variables
model.x = pyo.Var(within=pyo.NonNegativeReals)
model.y = pyo.Var(within=pyo.NonNegativeReals)

# 3. Define Objective (Minimize x + 2*y)
model.obj = pyo.Objective(expr=model.x + 2*model.y, sense=pyo.minimize)

# 4. Define Constraints
model.con1 = pyo.Constraint(expr=3*model.x + 4*model.y >= 12)
model.con2 = pyo.Constraint(expr=2*model.x + 5*model.y >= 10)

# 5. Solve
solver = pyo.SolverFactory('glpk')
results = solver.solve(model)

# 6. Access Results
print(f"x = {pyo.value(model.x)}, y = {pyo.value(model.y)}")
```

## Critical Rules

### ✅ DO

- **Scale your problem** - Like all optimization, Pyomo works best if variables and constraints are within a similar order of magnitude (e.g., 0.1 to 100).
- **Use pyo.value()** - Always wrap variables in `pyo.value(model.x)` to get their numerical result after solving.
- **Check Solver Status** - Verify the solver reached an optimal solution using `results.solver.termination_condition`.
- **Use pyo.Set and pyo.Param** - For complex models, organize data using these components rather than raw Python lists/dicts.
- **Specify Solver Options** - Use `solver.options['tm_limit'] = 60` (or similar depending on the solver) to manage execution time.
- **Use model.dual** - If the solver supports it, enable dual suffix to get shadow prices (important for economic analysis).

### ❌ DON'T

- **Don't use Python math functions** - Use `pyo.log()`, `pyo.exp()`, `pyo.sqrt()` instead of `math.log()` or `np.log()` inside expressions.
- **Don't use if statements in constraints** - Constraints must be algebraic. Use "Big-M" notation or logical reformulations for conditional logic.
- **Avoid "Hard" non-linearities** - Functions like `abs(x)` or `max(x, y)` are non-smooth and can break many solvers. Use reformulations.
- **Don't forget the Intercept/Constant** - Constraints like `x + y = 10` are fine, but ensure your units are consistent.

## Anti-Patterns (NEVER)

```python
import pyomo.environ as pyo
import numpy as np

# ❌ BAD: Using NumPy/Math functions in constraints
# model.con = pyo.Constraint(expr=np.sin(model.x) <= 0.5)

# ✅ GOOD: Use Pyomo-compatible functions
model.con = pyo.Constraint(expr=pyo.sin(model.x) <= 0.5)

# ❌ BAD: Using Python IF for conditional constraints
# if model.x > 10:
#     model.con = pyo.Constraint(expr=model.y <= 5)

# ✅ GOOD: Using Big-M notation (for binary variable z)
# y <= 5 + M * (1 - z)
# x >= 10 - M * (1 - z)

# ❌ BAD: Printing a variable directly
# print(model.x) # Returns a reference object, not a number!

# ✅ GOOD: Use value()
print(pyo.value(model.x))
```

## Modeling with Sets and Indices

### Large-Scale Model Construction

```python
model = pyo.ConcreteModel()

# Data
products = ['A', 'B', 'C']
profit = {'A': 10, 'B': 20, 'C': 15}
limit = 100

# Components
model.P = pyo.Set(initialize=products)
model.x = pyo.Var(model.P, within=pyo.NonNegativeReals)

# Indexed Objective
def obj_rule(model):
    return sum(profit[p] * model.x[p] for p in model.P)
model.obj = pyo.Objective(rule=obj_rule, sense=pyo.maximize)

# Indexed Constraint
def limit_rule(model):
    return sum(model.x[p] for p in model.P) <= limit
model.con = pyo.Constraint(rule=limit_rule)
```

## Non-Linear Programming (NLP)

### Using the IPOPT Solver

```python
# Solving: Minimize (x-2)^2 + (y-2)^2
model = pyo.ConcreteModel()
model.x = pyo.Var(initialize=0) # Initialization is CRUCIAL for NLP
model.y = pyo.Var(initialize=0)

model.obj = pyo.Objective(expr=(model.x - 2)**2 + (model.y - 2)**2)

# Constraint: x^2 + y <= 1
model.con = pyo.Constraint(expr=model.x**2 + model.y <= 1)

# Solve with IPOPT
solver = pyo.SolverFactory('ipopt')
solver.solve(model)
```

## Mixed-Integer Linear Programming (MILP)

### Discrete Decisions

```python
# Binary variable: 1 if we open a warehouse, 0 otherwise
model.use_warehouse = pyo.Var(within=pyo.Binary)

# Integer variable: Number of trucks to buy
model.num_trucks = pyo.Var(within=pyo.NonNegativeIntegers)

# Conditional logic: If warehouse is not used, trucks must be 0
# trucks <= Capacity * use_warehouse
model.cap_con = pyo.Constraint(expr=model.num_trucks <= 100 * model.use_warehouse)
```

## Practical Workflows

### 1. Diet Problem (Classic Linear Programming)

```python
def solve_diet(foods, nutrients, costs, requirements):
    model = pyo.ConcreteModel()
    model.F = pyo.Set(initialize=foods)
    model.N = pyo.Set(initialize=nutrients)
    
    model.x = pyo.Var(model.F, within=pyo.NonNegativeReals)
    
    # Minimize cost
    model.obj = pyo.Objective(expr=sum(costs[f] * model.x[f] for f in model.F))
    
    # Meet nutrient requirements
    def nutrient_rule(model, n):
        return sum(nutrients[f][n] * model.x[f] for f in model.F) >= requirements[n]
    model.con = pyo.Constraint(model.N, rule=nutrient_rule)
    
    pyo.SolverFactory('glpk').solve(model)
    return {f: pyo.value(model.x[f]) for f in model.F}
```

### 2. Blending Problem (Chemical/Process Engineering)

```python
# Balancing component fractions in a mixture
# Note: Often becomes non-linear (NLP) if both flow and fraction are variables
def blend_optimization(inputs, target_purity):
    model = pyo.ConcreteModel()
    # ... model setup ...
    # con: sum(flow[i] * purity[i]) / sum(flow[i]) == target_purity
    # becomes: sum(flow[i] * purity[i]) == target_purity * sum(flow[i]) (Linearized)
```

## Performance Optimization

### Solver Choice

- **GLPK/CBC**: Good for free/open-source LP/MIP.
- **Gurobi/CPLEX**: Industrial standards (extremely fast for large MIP).
- **IPOPT**: Best for smooth Non-linear (NLP).

### Warm Starts

For iterative optimizations, use