pyomo

# 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 the previous solution as a starting point.

```python
# For NLP solvers like IPOPT
model.x.set_value(prev_x_value)
solver.solve(model)
```

## Common Pitfalls and Solutions

### Termination Condition

Always check why the solver stopped.

```python
from pyomo.opt import TerminationCondition

results = solver.solve(model)
if results.solver.termination_condition == TerminationCondition.optimal:
    print("Success")
elif results.solver.termination_condition == TerminationCondition.infeasible:
    print("Check your constraints!")
```

### Non-Convexity in NLP

If your NLP model has multiple local minima, IPOPT might get stuck.

```python
# ✅ Solution: 
# 1. Provide multiple different initial guesses (multistart).
# 2. Use a global solver like BARON or SCIP.
```

### Indexing with Variables

`model.x[model.y]` where `y` is a `Var` is illegal.

```python
# ✅ Solution: Use model.AddElement or binary variable reformulations.
```

## Best Practices

1. Always scale your problem - variables and constraints should be within similar orders of magnitude
2. Use `pyo.value()` to extract numerical results from variables after solving
3. Check solver termination conditions to verify optimality or diagnose issues
4. Use `pyo.Set` and `pyo.Param` for organizing data in complex models
5. Initialize variables for NLP problems - good starting points are crucial
6. Use Pyomo-compatible functions (`pyo.log()`, `pyo.exp()`, etc.) instead of NumPy/math functions in expressions
7. Reformulate conditional logic using Big-M notation or binary variables
8. Avoid non-smooth functions like `abs()` or `max()` - use reformulations
9. Enable dual suffix for shadow prices when economic analysis is needed
10. Set solver time limits and options appropriately for your problem size

Pyomo is the ultimate tool for turning high-level mathematical abstractions into solved business and scientific problems. Its ability to bridge the gap between algebraic modeling and high-performance solvers makes it the foundation of modern prescriptive analytics.