Mathematical Optimization: Business Introduction#

What is Mathematical Optimization?#

Mathematical optimization (also called mathematical programming) is the process of finding the best solution from a set of possible choices, subject to constraints. Unlike simulation that explores “what if” scenarios, or heuristics that provide reasonable solutions, optimization mathematically guarantees the best achievable outcome given your objectives and limitations.

Why Use Mathematical Optimization?#

When Traditional Approaches Fall Short#

Many business problems involve making decisions with:

• Multiple competing objectives: Minimize cost while maximizing quality
• Complex constraints: Limited resources, capacity restrictions, regulatory requirements
• Combinatorial complexity: Thousands or millions of possible solutions to evaluate

For these problems, trial-and-error or even sophisticated heuristics may miss the optimal solution entirely. Mathematical optimization systematically searches the solution space to find provably optimal answers.

Core Concepts#

Objective Function#

What you want to optimize:

• Minimize: Costs, time, risk, waste, energy consumption
• Maximize: Revenue, throughput, efficiency, customer satisfaction

Decision Variables#

What you can control:

• Quantities to produce
• Resources to allocate
• Schedules to set
• Routes to follow

Constraints#

What limits your options:

• Budget limitations
• Capacity restrictions
• Time windows
• Physical limitations
• Regulatory requirements
• Business rules

When to Use Each Approach#

Mathematical Optimization#

Best for: Problems with clear objectives, well-defined constraints, and where you need the provably best solution.

Examples:

• Resource allocation across projects with budget constraints
• Production planning to minimize cost while meeting demand
• Portfolio optimization balancing risk and return
• Workforce scheduling with labor rules and coverage requirements
• Supply chain network design minimizing total cost

Advantages:

• Guarantees optimality (for solvable problem types)
• Handles complex constraint interactions automatically
• Provides sensitivity analysis (what happens if constraints change?)
• Scales to large problems with modern solvers

Limitations:

• Requires mathematically expressible objectives and constraints
• Some problem types (nonlinear, non-convex) are computationally hard
• May not capture all real-world complexities

Heuristics & Rules-Based Approaches#

Best for: Problems where optimization is too slow, or where simple rules capture sufficient domain knowledge.

Examples:

• First-come-first-served scheduling
• Greedy allocation algorithms
• Priority-based dispatching

Advantages:

• Fast to compute
• Easy to explain to stakeholders
• Can incorporate hard-to-quantify expertise

Limitations:

• No guarantee of quality
• May be far from optimal
• Difficult to handle constraint interactions

Simulation#

Best for: Understanding system behavior, quantifying uncertainty, or when the system is too complex to optimize directly.

Examples:

• Modeling customer flow through a service center
• Evaluating risk across uncertain scenarios
• Testing what-if scenarios

Advantages:

• Models complex, stochastic systems
• Evaluates variability and risk
• No need for closed-form equations

Limitations:

• Doesn’t find optimal decisions (only evaluates given scenarios)
• Can be computationally expensive
• Requires many simulation runs for statistical significance

Simulation + Optimization#

Best for: Complex systems where the objective can’t be written analytically, but can be evaluated through simulation.

Examples:

• Optimizing inventory policies in a supply chain with complex lead times
• Finding best control parameters for a manufacturing process
• Tuning staffing levels in a service system

Approach: Use simulation to evaluate objective function, use optimization to search for best parameters.

Generic Business Value#

Mathematical optimization delivers value across industries:

Cost Reduction#

• Minimize production costs while meeting demand
• Reduce transportation and logistics expenses
• Optimize energy consumption
• Minimize inventory holding costs

Revenue Maximization#

• Optimal pricing strategies
• Product mix optimization
• Capacity allocation to high-value customers
• Resource allocation to high-return opportunities

Resource Efficiency#

• Maximize throughput from limited resources
• Optimal workforce scheduling and assignment
• Equipment utilization optimization
• Space utilization in facilities

Service Quality#

• Minimize customer wait times
• Balance workload across resources
• Meet service level agreements with minimal cost
• Optimize coverage and availability

Risk Management#

• Portfolio diversification
• Robust solutions under uncertainty
• Minimize worst-case scenarios
• Balance risk and return

Generic Examples Across Industries#

Resource Allocation#

Problem: Allocate limited budget across N projects to maximize total return, where projects have different costs and returns, and some projects have dependencies.

Optimization approach: Mixed-integer programming to select project portfolio maximizing total value subject to budget and dependency constraints.

Scheduling with Precedence#

Problem: Schedule tasks with precedence constraints (some tasks must finish before others can start) to minimize total completion time.

Optimization approach: Constraint programming or mixed-integer programming with precedence and resource capacity constraints.

Portfolio Optimization#

Problem: Select mix of assets that maximizes expected return while keeping risk below a threshold, with constraints on diversification.

Optimization approach: Quadratic programming (minimize variance of portfolio return) or robust optimization (worst-case scenarios).

Network Flow Optimization#

Problem: Route flow through a network (supply chain, communication network, transportation) to minimize cost while satisfying demand at destinations.

Optimization approach: Linear programming (network flow problems have special structure enabling efficient solution).

Cutting Stock / Bin Packing#

Problem: Cut raw materials into required pieces minimizing waste, or pack items into containers minimizing number of containers used.

Optimization approach: Mixed-integer programming with complex cutting patterns or packing configurations.

Modern Tooling Landscape#

Python has become the dominant platform for optimization due to:

• Rich ecosystem: Multiple modeling languages and solver interfaces
• Integration: Easy connection to data pipelines, machine learning, simulation
• Accessibility: Free, open-source tools competitive with commercial software
• Flexibility: From rapid prototyping to production deployment

The optimization landscape spans:

• Open-source solvers: Free, high-quality solvers (HiGHS, SCIP, IPOPT)
• Algebraic modeling: Express models in mathematical notation (Pyomo, CVXPY, PuLP)
• Commercial solvers: Premium performance for large-scale problems (Gurobi, CPLEX)
• Specialized libraries: Multi-objective optimization, dynamic optimization, constraint programming

Getting Started#

1. Identify your problem type: Linear? Integer variables? Nonlinear? Convex?
2. Start simple: Begin with small models, validate results
3. Iterate: Add complexity gradually, test constraint impacts
4. Validate: Compare optimization results to current practice or heuristics
5. Sensitivity analysis: Understand how solution changes with parameters

Mathematical optimization is a powerful tool when applied to well-defined problems. The key is understanding when optimization is the right approach, what problem type you’re solving, and selecting appropriate tools for your use case.

Rapid Discovery: What is Mathematical Optimization?#

Technical Overview#

Mathematical optimization (mathematical programming) is the discipline of finding the best element from a feasible set with respect to an objective criterion. Formally:

minimize (or maximize)    f(x)
subject to                g_i(x) ≤ 0,  i = 1,...,m
                          h_j(x) = 0,  j = 1,...,p
                          x ∈ X

Where:

• f(x): Objective function to minimize/maximize
• x: Decision variables (what we control)
• g_i(x): Inequality constraints
• h_j(x): Equality constraints
• X: Domain of variables (e.g., real numbers, integers)

Problem Type Taxonomy#

The structure of f(x), g_i(x), h_j(x), and X determines problem type, which drives solver selection.

Linear Programming (LP)#

Definition: Linear objective, linear constraints, continuous variables.

minimize    c^T x
subject to  Ax ≤ b
            x ≥ 0

Example: Maximize profit from production of multiple products with resource constraints.

Characteristics:

• Polynomial-time solvable (simplex, interior point methods)
• Global optimum guaranteed
• Well-understood theory
• Efficient solvers (HiGHS, GLPK, commercial)

When to use: Resource allocation, blending problems, network flows, diet problems, transportation problems.

Mixed-Integer Linear Programming (MILP/MIP)#

Definition: LP where some/all variables must be integers.

minimize    c^T x
subject to  Ax ≤ b
            x_i ∈ {0,1} or x_i ∈ ℤ for some i

Example: Select which facilities to open (binary decisions) to minimize total cost.

Characteristics:

• NP-hard (computationally difficult)
• Branch-and-bound, branch-and-cut algorithms
• Modern solvers extremely effective for practical instances
• Special structure (e.g., network flows) helps

When to use: Yes/no decisions (selection, assignment), scheduling with discrete time periods, facility location, cutting stock, bin packing.

Quadratic Programming (QP)#

Definition: Quadratic objective, linear constraints.

minimize    (1/2) x^T Q x + c^T x
subject to  Ax ≤ b

Example: Portfolio optimization minimizing variance (quadratic) subject to return requirements.

Characteristics:

• If Q is positive semidefinite (convex), efficiently solvable
• If Q is indefinite (non-convex), NP-hard
• Interior point methods effective for convex case

When to use: Portfolio optimization, least-squares problems, model predictive control.

Nonlinear Programming (NLP)#

Definition: Nonlinear objective or constraints.

minimize    f(x)
subject to  g_i(x) ≤ 0

Example: Minimize chemical reactor cost (nonlinear engineering equations) subject to safety constraints.

Characteristics:

• Wide range of difficulty depending on problem structure
• Local optimum vs global optimum distinction critical
• Gradient-based methods (BFGS, SQP, interior point)
• Derivative-free methods for non-smooth functions

When to use: Engineering design, parameter estimation, fitting nonlinear models, physical system optimization.

Convex Optimization#

Definition: Convex objective function, convex feasible region.

A function f is convex if:

f(αx + (1-α)y) ≤ αf(x) + (1-α)f(y)  for all α ∈ [0,1]

Example: Minimize convex loss function subject to norm constraints.

Characteristics:

• Any local optimum is global optimum
• Polynomial-time solvable for many subclasses
• LP and convex QP are special cases
• Includes conic programming (SOCP, SDP)

When to use: Machine learning (SVM, logistic regression), signal processing, control systems, robust optimization.

Constraint Programming (CP)#

Definition: Logical and combinatorial constraints, often with discrete variables.

Example: Scheduling tasks with precedence constraints, resource limits, and no-overlap requirements.

Characteristics:

• Domain propagation and constraint inference
• Effective for scheduling, timetabling, configuration
• Different paradigm from mathematical programming

When to use: Scheduling with complex logical rules, rostering, configuration problems, puzzles.

Multi-Objective Optimization#

Definition: Multiple objective functions simultaneously.

minimize    [f_1(x), f_2(x), ..., f_k(x)]
subject to  g_i(x) ≤ 0

Example: Minimize cost AND minimize emissions (competing objectives).

Characteristics:

• No single “optimal” solution
• Find Pareto frontier (solutions where improving one objective worsens another)
• Requires preference specification or generate full frontier

When to use: Trade-off analysis, multi-criteria decision making, design optimization with competing goals.

Solver vs Modeling Language Distinction#

Solver#

The engine that implements optimization algorithms and computes solutions.

Examples:

• Open-source: HiGHS, CBC, GLPK, SCIP, IPOPT, OSQP
• Commercial: Gurobi, CPLEX, MOSEK, XPRESS

Role: Takes problem in standard form (MPS, LP file, or direct API) and returns solution.

Modeling Language#

The interface for expressing optimization problems in human-readable form.

Examples:

• Python: Pyomo, CVXPY, PuLP, OR-Tools
• Others: AMPL, GAMS, JuMP (Julia)

Role: Translates user model to solver input format, manages data, retrieves results.

Analogy: Modeling language is like SQL (express what you want), solver is like database engine (compute the result).

Optimization Algorithm Approaches#

Gradient-Based Methods#

Principle: Use derivatives to determine descent direction.

Algorithms:

• Gradient descent: Move in direction of negative gradient
• Newton’s method: Use second derivatives (Hessian) for faster convergence
• Quasi-Newton (BFGS): Approximate Hessian efficiently
• Sequential Quadratic Programming (SQP): Solve sequence of QP subproblems
• Interior point methods: Move through interior of feasible region

Requirements: Objective and constraints must be differentiable.

Advantages:

• Fast convergence for smooth problems
• Efficient for large-scale problems
• Mature theory and implementations

Limitations:

• May converge to local optima (non-convex problems)
• Requires gradient computation (analytical or automatic differentiation)

When to use: Continuous optimization, smooth functions, large scale.

Derivative-Free Methods#

Principle: Optimize without computing gradients.

Algorithms:

• Nelder-Mead simplex: Geometric search using simplex shape
• Powell’s method: Sequential line searches
• Pattern search: Evaluate on geometric patterns
• Trust region methods: Sample within trusted radius

Advantages:

• Handle non-smooth, noisy objectives
• Simple to implement
• No derivative computation needed

Limitations:

• Slower convergence than gradient-based
• Struggle with high dimensions
• No guarantee of global optimum

When to use: Black-box objectives (simulation-based), non-smooth functions, small to medium dimension.

Combinatorial/Discrete Optimization#

Principle: Systematically search discrete solution space.

Algorithms:

• Branch-and-bound: Partition space, bound subproblem solutions
• Branch-and-cut: Branch-and-bound with cutting planes
• Dynamic programming: Solve subproblems recursively
• Constraint propagation: Eliminate infeasible values

Characteristics:

• Exact methods guarantee optimal solution
• Computational complexity often exponential
• Modern implementations highly optimized

When to use: Integer variables, combinatorial structure.

Metaheuristics#

Principle: Heuristic search strategies inspired by natural processes.

Algorithms:

• Genetic algorithms: Evolutionary selection and crossover
• Simulated annealing: Probabilistic acceptance of worse solutions
• Particle swarm optimization: Population-based search
• Tabu search: Guided local search with memory

Advantages:

• Applicable to any problem type
• Can escape local optima
• Handle large, complex, black-box problems

Limitations:

• No optimality guarantee
• Many parameters to tune
• Computational cost (many evaluations)

When to use: Complex, non-convex problems where exact methods too slow, multi-objective optimization, robust optimization.

Key Insights for Solver Selection#

1. Problem type drives solver choice: LP solvers can’t handle integers, MILP solvers can’t handle nonlinearity. Know your problem type.

2. Convexity is gold: Convex problems have global optima and efficient algorithms. Test if your problem is convex.

3. Size matters differently: LP scales to millions of variables. MILP scales to thousands to hundreds of thousands. NLP scaling depends heavily on structure.

4. Modeling language adds convenience: For rapid development, algebraic modeling (Pyomo, CVXPY) beats direct API coding.

5. Commercial vs open-source gap narrowing: Open-source solvers (HiGHS, SCIP) competitive for many problems. Commercial edge on largest MILP instances.

6. Solver backend matters more than frontend: For production performance, solver (Gurobi vs CBC) dominates modeling language choice (Pyomo vs PuLP).

Problem Type Decision Tree#

Is objective and constraints linear?
├─ YES, all continuous → LP
│  └─ Solvers: HiGHS, GLPK, Gurobi, CPLEX
├─ YES, some integer → MILP
│  └─ Solvers: CBC, SCIP, HiGHS, Gurobi, CPLEX
└─ NO (nonlinear)
   ├─ Is problem convex?
   │  ├─ YES → Convex optimization
   │  │  └─ Solvers: ECOS, SCS, MOSEK, Clarabel
   │  └─ NO → General NLP
   │     └─ Solvers: IPOPT, SNOPT, Knitro
   └─ Logical constraints, scheduling? → Constraint Programming
      └─ Solvers: OR-Tools CP-SAT, Gecode

Multiple objectives?
└─ Multi-objective optimization (pymoo, weighted sum, epsilon-constraint)

Dynamic over time with differential equations?
└─ Dynamic optimization (GEKKO, Pyomo.DAE)

Next Steps in S1-Rapid#

The following pages investigate specific Python libraries:

• scipy.optimize: Built-in Python scientific computing
• OR-Tools: Google’s operations research toolkit
• Pyomo: Algebraic modeling language
• CVXPY: Convex optimization specialist
• PuLP: Simple LP/MILP modeling

Each library page documents:

• Problem types supported
• Installation and basic example
• API approach and design philosophy
• Solver backends available
• Maturity and community indicators
• When to choose this library

CVXPY: Convex Optimization in Python#

Overview#

CVXPY is a Python-embedded modeling language for convex optimization problems. It provides automatic verification of problem convexity through Disciplined Convex Programming (DCP) analysis, ensuring that your problem can be solved globally and efficiently.

Key finding: CVXPY is the specialist for convex optimization, with Stanford academic origins and JMLR publication (2016). Unlike general-purpose tools, CVXPY verifies convexity at modeling time, preventing common errors.

Installation#

pip install cvxpy

CVXPY includes open-source solvers (Clarabel, SCS, OSQP, ECOS) by default. Additional solvers (MOSEK, Gurobi, CPLEX) can be used if installed.

Problem Types Supported#

CVXPY is specialized for convex optimization:

Core insight: If your problem is convex, CVXPY is likely the best choice. DCP analysis catches modeling errors early.

What is Convex Optimization?#

A problem is convex if:

1. Objective function is convex (for minimization)
2. Feasible region is convex

Why convexity matters:

• Any local optimum is a global optimum
• Polynomial-time algorithms available
• No initialization sensitivity
• Reliable, predictable solving

Common convex problems:

• Linear programming
• Least-squares regression
• Lasso, ridge regression
• Support vector machines (SVM)
• Portfolio optimization (mean-variance)
• Signal processing (filter design)
• Control systems (LQR, MPC)

Disciplined Convex Programming (DCP)#

DCP is a ruleset for constructing convex optimization problems. CVXPY automatically checks DCP compliance.

DCP rules:

• Convex functions can only be minimized or bounded above
• Concave functions can only be maximized or bounded below
• Affine functions can be used anywhere
• Composition rules ensure convexity is preserved

Example of DCP violation:

import cvxpy as cp

x = cp.Variable()
# ERROR: Cannot minimize a concave function
objective = cp.Minimize(-cp.sqrt(x))  # sqrt is concave

CVXPY raises an error immediately, preventing silent failures.

Basic Example: Linear Programming#

import cvxpy as cp

# Variables
x = cp.Variable()
y = cp.Variable()

# Objective: minimize -x - 2y (equivalent to maximize x + 2y)
objective = cp.Minimize(-x - 2*y)

# Constraints
constraints = [
    2*x + y <= 20,
    x + 2*y <= 20,
    x >= 0,
    y >= 0
]

# Create and solve problem
problem = cp.Problem(objective, constraints)
problem.solve(solver=cp.CLARABEL)

print(f"Status: {problem.status}")
print(f"Optimal value: {-problem.value}")  # Negate for maximization
print(f"x = {x.value}")
print(f"y = {y.value}")

Basic Example: Quadratic Programming (Portfolio Optimization)#

import cvxpy as cp
import numpy as np

# Portfolio optimization: minimize risk subject to return target
n_assets = 4
returns = np.array([0.10, 0.15, 0.12, 0.08])  # Expected returns
cov_matrix = np.array([  # Covariance matrix (risk)
    [0.04, 0.01, 0.02, 0.01],
    [0.01, 0.06, 0.01, 0.02],
    [0.02, 0.01, 0.05, 0.01],
    [0.01, 0.02, 0.01, 0.03]
])

# Variables: portfolio weights
w = cp.Variable(n_assets)

# Objective: minimize portfolio variance (risk)
risk = cp.quad_form(w, cov_matrix)
objective = cp.Minimize(risk)

# Constraints
constraints = [
    cp.sum(w) == 1,           # Weights sum to 1
    w >= 0,                    # No short selling
    returns @ w >= 0.11        # Target return of 11%
]

problem = cp.Problem(objective, constraints)
problem.solve()

print(f"Optimal portfolio: {w.value}")
print(f"Expected return: {returns @ w.value:.2%}")
print(f"Portfolio variance: {problem.value:.4f}")
print(f"Portfolio std dev: {np.sqrt(problem.value):.2%}")

Basic Example: L1 Regression (Lasso)#

import cvxpy as cp
import numpy as np

# Generate synthetic data
np.random.seed(1)
n, p = 100, 20
X = np.random.randn(n, p)
beta_true = np.zeros(p)
beta_true[:5] = np.array([1, -2, 0.5, -1, 3])  # Sparse true coefficients
y = X @ beta_true + 0.1 * np.random.randn(n)

# Lasso regression: minimize ||Xβ - y||² + λ||β||₁
beta = cp.Variable(p)
lambda_reg = 0.1

objective = cp.Minimize(cp.sum_squares(X @ beta - y) + lambda_reg * cp.norm(beta, 1))
problem = cp.Problem(objective)
problem.solve()

print(f"True coefficients: {beta_true[:5]}")
print(f"Estimated coefficients: {beta.value[:5]}")
print(f"Number of non-zero: {np.sum(np.abs(beta.value) > 0.01)}")

Basic Example: Constraint with Norms (SOCP)#

import cvxpy as cp

# Minimize ||x - center||₂ subject to linear constraints
x = cp.Variable(3)
center = np.array([1, 2, 3])

objective = cp.Minimize(cp.norm(x - center, 2))

constraints = [
    cp.sum(x) == 9,
    x >= 0
]

problem = cp.Problem(objective, constraints)
problem.solve()

print(f"Closest point: {x.value}")
print(f"Distance: {problem.value}")

API Design Philosophy#

Algebraic modeling with DCP verification: Express problems naturally while CVXPY ensures convexity.

Advantages:

• Automatic convexity verification: Catch modeling errors at construction time
• Readable code: Mathematical notation
• Solver abstraction: CVXPY reformulates problem for solver capabilities
• Parameters: Separate data from structure, enable rapid re-solves

Disadvantages:

• Limited to convex problems: Cannot handle general nonlinear
• No integer variables (or very limited support via MI-CVXPY)
• DCP learning curve: Must understand convexity rules

Solver Backends#

CVXPY automatically selects appropriate solver based on problem type.

Bundled (Free):#

• Clarabel: Rust-based interior point (default for many problems)
• SCS: Splitting Conic Solver (robust, medium accuracy)
• OSQP: Operator Splitting QP (fast for QP)
• ECOS: Embedded Conic Solver (small to medium SOCP)
• CVXOPT: Python convex optimization

Commercial (if installed):#

• MOSEK: High-performance conic (SOCP, SDP)
• Gurobi: QP and conic support
• CPLEX: QP support

Specialized:#

• GLPK, CBC: LP/MILP via Pyomo interface
• NAG: Numerical Algorithms Group

Specify solver with problem.solve(solver=cp.MOSEK).

When to Use CVXPY#

✅ Good fit when:#

• Problem is convex: LP, convex QP, SOCP, SDP, convex constraints
• Want convexity verification: Prevent modeling errors with DCP
• Machine learning: Lasso, ridge, SVM, logistic regression
• Portfolio optimization: Mean-variance, risk parity
• Signal processing: Filter design, compressed sensing
• Control theory: LQR, MPC with convex constraints
• Parameter tuning: Rapid re-solves with different data (cp.Parameter)

❌ Not suitable when:#

• Integer variables: No MILP support → use PuLP, Pyomo, OR-Tools
• Nonconvex NLP: Non-convex objectives/constraints → use scipy, Pyomo+IPOPT
• Combinatorial optimization: Scheduling, routing → use OR-Tools
• Large-scale MILP: Commercial solvers faster → Gurobi, CPLEX directly

Community and Maturity#

CVXPY is the standard tool for convex optimization in Python, widely used in academia and industry.

Key Findings from Research#

1. DCP verification is killer feature: Automatically checks convexity at modeling time. No other Python library does this.

2. Academic pedigree: Developed by Stephen Boyd’s group at Stanford, based on CVX (MATLAB). Used in teaching convex optimization.

3. Automatic reformulation: CVXPY transforms your problem into standard conic form for solvers. You model naturally, CVXPY handles the details.

4. Parameters enable fast re-solves: Declare cp.Parameter for data that changes, dramatically faster than rebuilding model.

5. Disciplined Convex Programming (DCP): Ruleset ensures convexity. Learning curve, but catches errors early.

Comparison with Alternatives#

Example Use Cases (Generic)#

Resource Allocation with Quadratic Cost#

Allocate resources to minimize quadratic cost (diminishing returns).

import cvxpy as cp

# Allocate budget across N options with diminishing returns
x = cp.Variable(n_options)

# Quadratic cost (diminishing returns modeled as x^2)
objective = cp.Minimize(cp.sum_squares(x))

constraints = [
    cp.sum(x) == total_budget,
    x >= min_allocation,
    x <= max_allocation
]

problem = cp.Problem(objective, constraints)
problem.solve()

Robust Optimization#

Minimize worst-case cost under uncertainty.

import cvxpy as cp

# Worst-case cost over uncertain parameters
x = cp.Variable(n)
u = cp.Variable(m)  # Uncertain parameters

# Minimize worst-case objective
objective = cp.Minimize(cp.max(cost_scenarios @ x))

constraints = [
    A @ x + B @ u <= b,
    cp.norm(u, 2) <= uncertainty_radius,  # Uncertainty set
    x >= 0
]

problem = cp.Problem(objective, constraints)
problem.solve()

Fair Resource Allocation#

Maximize minimum allocation (max-min fairness).

import cvxpy as cp

# Allocate resources fairly (maximize minimum share)
x = cp.Variable(n_users)
min_share = cp.Variable()

objective = cp.Maximize(min_share)

constraints = [
    x >= min_share,  # Each user gets at least min_share
    cp.sum(x) <= total_capacity,
    x >= 0
]

problem = cp.Problem(objective, constraints)
problem.solve()

Advanced Features#

Parameters for Fast Re-solving#

# Define parameter (data that changes)
lambda_param = cp.Parameter(nonneg=True)
beta = cp.Variable(p)

objective = cp.Minimize(cp.sum_squares(X @ beta - y) + lambda_param * cp.norm(beta, 1))
problem = cp.Problem(objective)

# Solve for different lambda values (fast!)
for lam in [0.01, 0.1, 1.0]:
    lambda_param.value = lam
    problem.solve(warm_start=True)
    print(f"Lambda={lam}: objective={problem.value}")

Constraints as Objects#

# Store constraint for later modification
budget_constraint = cp.sum(x) <= budget_param
problem = cp.Problem(objective, [budget_constraint, ...])

# Later, change budget and re-solve
budget_param.value = new_budget
problem.solve()

Dual Variables#

# Access dual variables (shadow prices)
problem.solve()
print(f"Dual value: {constraints[0].dual_value}")

Problem Classes CVXPY Recognizes#

CVXPY automatically classifies problems for solver selection:

• LP: Linear Programming
• QP: Quadratic Programming
• SOCP: Second-Order Cone Programming (norm constraints)
• SDP: Semidefinite Programming (matrix constraints)
• ExpCone: Exponential cone (entropy, log)
• PowCone: Power cone (geometric programming)

Check with: print(problem.classification)

References#

• Official documentation: https://www.cvxpy.org
• GitHub: https://github.com/cvxpy/cvxpy
• Paper: Diamond & Boyd, “CVXPY: A Python-Embedded Modeling Language for Convex Optimization”, JMLR 2016
• Tutorial: https://www.cvxpy.org/tutorial/
• Book: Boyd & Vandenberghe, “Convex Optimization” (free online)

Summary#

CVXPY is the specialist for convex optimization:

• Automatic convexity verification (DCP analysis)
• Natural mathematical notation
• Excellent for machine learning, portfolio optimization, signal processing
• Stanford academic origins, widely used in education and research

Choose CVXPY when:

• Your problem is convex (or you want to verify it is)
• Machine learning applications (Lasso, SVM, logistic regression)
• Portfolio optimization, control theory, signal processing
• You want to catch modeling errors early with DCP

Look elsewhere for:

• Integer variables (PuLP, Pyomo, OR-Tools)
• General nonlinear (scipy, Pyomo)
• Large-scale MILP (Gurobi, CPLEX)

If your problem is convex, start with CVXPY. The DCP verification alone is worth it.

Google OR-Tools: Industrial-Strength Operations Research#

Overview#

OR-Tools is Google’s open-source operations research toolkit. It provides high-performance solvers for constraint programming (CP), linear/integer programming (LP/MILP), vehicle routing, and graph algorithms. Written in C++ with Python, Java, and C# bindings.

Key finding: With 12.6k GitHub stars, OR-Tools is the most popular open-source optimization framework. Actively maintained by Google with Python 3.13 support as of 2025.

Installation#

pip install ortools

OR-Tools is self-contained with solvers bundled (no external solver installation required).

Problem Types Supported#

Basic Example: Linear Programming#

from ortools.linear_solver import pywraplp

# Create solver instance
solver = pywraplp.Solver.CreateSolver('GLOP')  # Google's LP solver

# Create variables: 0 <= x <= infinity, 0 <= y <= infinity
x = solver.NumVar(0, solver.infinity(), 'x')
y = solver.NumVar(0, solver.infinity(), 'y')

# Create constraints
# Constraint 1: 2x + y <= 20
solver.Add(2 * x + y <= 20)

# Constraint 2: x + 2y <= 20
solver.Add(x + 2 * y <= 20)

# Objective: maximize x + 2y
solver.Maximize(x + 2 * y)

# Solve
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print(f'Optimal solution found:')
    print(f'x = {x.solution_value()}')
    print(f'y = {y.solution_value()}')
    print(f'Optimal value = {solver.Objective().Value()}')
else:
    print('No optimal solution found.')

Basic Example: Mixed-Integer Programming#

from ortools.linear_solver import pywraplp

# Create MILP solver (using SCIP backend)
solver = pywraplp.Solver.CreateSolver('SCIP')

# Integer variables: x, y ∈ {0, 1, 2, ...}
x = solver.IntVar(0, solver.infinity(), 'x')
y = solver.IntVar(0, solver.infinity(), 'y')

# Constraints
solver.Add(2 * x + y <= 20)
solver.Add(x + 2 * y <= 20)

# Objective: maximize x + 2y
solver.Maximize(x + 2 * y)

# Solve
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL:
    print(f'x = {x.solution_value()}')  # Integer values
    print(f'y = {y.solution_value()}')
    print(f'Optimal value = {solver.Objective().Value()}')

Basic Example: Constraint Programming (CP-SAT)#

CP-SAT is Google’s state-of-the-art constraint programming solver, excellent for scheduling and combinatorial problems.

from ortools.sat.python import cp_model

# Create model
model = cp_model.CpModel()

# Variables with domain [0, 10]
x = model.NewIntVar(0, 10, 'x')
y = model.NewIntVar(0, 10, 'y')
z = model.NewIntVar(0, 10, 'z')

# Constraint: x != y (all different)
model.Add(x != y)
model.Add(y != z)
model.Add(x != z)

# Constraint: x + y + z == 15
model.Add(x + y + z == 15)

# Objective: maximize x * y * z
model.Maximize(x * y * z)

# Solve
solver = cp_model.CpSolver()
status = solver.Solve(model)

if status == cp_model.OPTIMAL:
    print(f'x = {solver.Value(x)}')
    print(f'y = {solver.Value(y)}')
    print(f'z = {solver.Value(z)}')
    print(f'Objective = {solver.ObjectiveValue()}')

Modules and Components#

OR-Tools is organized into specialized modules:

1. Linear/Integer Programming (linear_solver)#

• GLOP: Google’s LP solver (primal/dual simplex)
• SCIP, CBC, GLPK: Bundled MILP solvers
• External: Can interface with Gurobi, CPLEX if installed

2. Constraint Programming (sat)#

• CP-SAT: Modern SAT-based CP solver
• Excellent for scheduling, assignment, sequencing
• Exploits Boolean satisfiability techniques
• Winner of multiple MiniZinc Challenge medals

3. Vehicle Routing (routing)#

• Specialized solvers for VRP variants
• Rich constraint modeling (time windows, capacities, multiple vehicles)
• Metaheuristics tuned for routing

4. Graph Algorithms (graph)#

• Min-cost flow, max-flow
• Shortest paths
• Assignment problems
• Network optimization

API Design Philosophy#

Object-oriented, solver-specific APIs: Different modules have different APIs (linear_solver vs CP-SAT vs routing).

Advantages:

• High performance (thin Python wrapper over C++)
• Rich features specific to problem type
• Bundled solvers (easy installation)

Disadvantages:

• Different APIs for different problem types (less unified than Pyomo)
• More verbose than algebraic modeling
• Harder to switch between solver types

When to Use OR-Tools#

✅ Good fit when:#

• Constraint programming: Scheduling, sequencing, combinatorial problems → CP-SAT
• Vehicle routing: Delivery routing, technician dispatch → routing module
• Production environment: Stable, high-performance, Google-backed
• Avoid external dependencies: Bundled solvers, no license management
• Integer programming: Strong MILP capability via SCIP
• Graph/network problems: Specialized efficient algorithms

❌ Not suitable when:#

• Nonlinear programming: No NLP support → use scipy, Pyomo+IPOPT, GEKKO
• Convex optimization: No specialized convex solvers → use CVXPY
• Algebraic modeling preference: More verbose than Pyomo/CVXPY
• Academic flexibility: Pyomo supports more solver backends

Solver Backends Available#

Bundled with OR-Tools:#

• GLOP: Google’s LP solver
• CP-SAT: Google’s CP solver
• SCIP: Open-source MILP (Apache 2.0 since 9.0)
• GLPK: GNU Linear Programming Kit
• CBC: COIN-OR Branch and Cut

External (if installed):#

• Gurobi: Commercial MILP
• CPLEX: IBM commercial MILP
• XPRESS: FICO commercial MILP

Use CreateSolver('SOLVER_NAME') to select backend.

Community and Maturity#

OR-Tools is production-grade and used extensively at Google and beyond.

Key Findings from Research#

1. CP-SAT dominance: Google’s CP-SAT solver wins MiniZinc Challenge medals, competitive with commercial CP solvers. Best-in-class for scheduling.

2. Self-contained: All solvers bundled. No license management, no external compilation. Major usability advantage.

3. Multiple APIs: Different modules (linear_solver, sat, routing) have distinct APIs. Less unified than Pyomo but more specialized.

4. Python 3.13 support: As of 2025, OR-Tools supports latest Python versions, showing active maintenance.

5. Vehicle routing specialization: Dedicated routing library with domain-specific constraints (time windows, capacities, precedence). Unique among Python optimization libraries.

Comparison with Alternatives#

Example Use Cases (Generic)#

Resource Assignment with Constraints#

Assign tasks to workers subject to skill requirements, availability, workload limits.

from ortools.sat.python import cp_model

model = cp_model.CpModel()

# Binary variables: worker_task[w, t] = 1 if worker w assigned to task t
worker_task = {}
for w in workers:
    for t in tasks:
        worker_task[w, t] = model.NewBoolVar(f'worker_{w}_task_{t}')

# Each task assigned to exactly one worker
for t in tasks:
    model.Add(sum(worker_task[w, t] for w in workers) == 1)

# Worker workload limits
for w in workers:
    model.Add(sum(duration[t] * worker_task[w, t] for t in tasks) <= max_workload[w])

# Minimize total cost
model.Minimize(sum(cost[w, t] * worker_task[w, t]
                   for w in workers for t in tasks))

solver = cp_model.CpSolver()
solver.Solve(model)

Scheduling with No-Overlap#

Schedule jobs on machines with no time overlap.

from ortools.sat.python import cp_model

model = cp_model.CpModel()

# Interval variables for each job
intervals = []
for j in jobs:
    start = model.NewIntVar(0, horizon, f'start_{j}')
    duration = durations[j]
    end = model.NewIntVar(0, horizon, f'end_{j}')
    interval = model.NewIntervalVar(start, duration, end, f'interval_{j}')
    intervals.append(interval)

# No overlap constraint
model.AddNoOverlap(intervals)

# Minimize makespan (max end time)
makespan = model.NewIntVar(0, horizon, 'makespan')
for j in jobs:
    model.Add(makespan >= intervals[j].EndExpr())

model.Minimize(makespan)

Network Flow Problem#

Route flow through network minimizing cost.

from ortools.graph import pywrapgraph

# Create min-cost flow solver
min_cost_flow = pywrapgraph.SimpleMinCostFlow()

# Add arcs: (source, destination, capacity, unit_cost)
for source, dest, capacity, cost in arcs:
    min_cost_flow.AddArcWithCapacityAndUnitCost(source, dest, capacity, cost)

# Set supply/demand at nodes
for node, supply in supplies.items():
    min_cost_flow.SetNodeSupply(node, supply)

# Solve
status = min_cost_flow.Solve()

if status == min_cost_flow.OPTIMAL:
    print(f'Minimum cost: {min_cost_flow.OptimalCost()}')
    for arc in range(min_cost_flow.NumArcs()):
        if min_cost_flow.Flow(arc) > 0:
            print(f'Arc {arc}: flow = {min_cost_flow.Flow(arc)}')

Performance Characteristics#

• LP (GLOP): Competitive with open-source solvers (GLPK, CLP), somewhat slower than commercial
• MILP (SCIP): Good performance, especially for structured problems
• CP-SAT: Excellent for scheduling and combinatorial problems, award-winning
• Routing: Highly optimized metaheuristics for vehicle routing variants

For maximum MILP performance on large instances, commercial solvers (Gurobi, CPLEX) still faster, but OR-Tools competitive for many practical problems.

References#

• Official documentation: https://developers.google.com/optimization
• GitHub repository: https://github.com/google/or-tools (12.6k stars)
• Examples: https://github.com/google/or-tools/tree/stable/examples/python
• CP-SAT: https://developers.google.com/optimization/cp/cp_solver

Summary#

OR-Tools is a comprehensive, production-ready optimization toolkit with:

• Best-in-class constraint programming (CP-SAT)
• Specialized vehicle routing solvers
• Bundled MILP solvers (SCIP, CBC, GLPK)
• Google backing and active development

Choose OR-Tools for:

• Scheduling, sequencing, combinatorial problems (CP-SAT)
• Vehicle routing and logistics
• Production environments requiring stable, maintained library
• Projects avoiding external solver installation

Look elsewhere for nonlinear programming, convex optimization, or algebraic modeling preference.

PuLP: Linear Programming in Python#

Overview#

PuLP is a linear and mixed-integer programming modeler written in Python. It is designed to be simple and accessible for LP/MILP problems, with an algebraic modeling interface and the CBC solver bundled by default.

Key finding: PuLP 3.3.0 (Sept 2025) offers the easiest entry point for LP/MILP optimization in Python. Part of COIN-OR project with MIT license. CBC solver included means zero external dependencies for basic use.

Installation#

pip install pulp

The CBC (COIN-OR Branch and Cut) solver is bundled, so you can solve MILP problems immediately without additional installation.

Problem Types Supported#

PuLP is focused exclusively on linear programming:

Philosophy: Do one thing well. PuLP targets LP/MILP with simple, clean API.

Basic Example: Linear Programming#

from pulp import *

# Create problem
prob = LpProblem("Example", LpMaximize)

# Variables
x = LpVariable("x", lowBound=0)
y = LpVariable("y", lowBound=0)

# Objective: maximize x + 2y
prob += x + 2*y, "Objective"

# Constraints
prob += 2*x + y <= 20, "Constraint1"
prob += x + 2*y <= 20, "Constraint2"

# Solve with bundled CBC solver
prob.solve()

# Results
print(f"Status: {LpStatus[prob.status]}")
print(f"Optimal value: {value(prob.objective)}")
print(f"x = {value(x)}")
print(f"y = {value(y)}")

Basic Example: Integer Programming#

from pulp import *

prob = LpProblem("Integer_Example", LpMaximize)

# Integer variables
x = LpVariable("x", lowBound=0, cat='Integer')
y = LpVariable("y", lowBound=0, cat='Integer')

# Objective
prob += x + 2*y

# Constraints
prob += 2*x + y <= 20
prob += x + 2*y <= 20

# Solve
prob.solve(PULP_CBC_CMD(msg=0))  # msg=0 suppresses solver output

print(f"x = {value(x)}")  # Integer value
print(f"y = {value(y)}")  # Integer value
print(f"Optimal value: {value(prob.objective)}")

Basic Example: Binary Selection Problem#

from pulp import *

# Select subset of projects to maximize value subject to budget
projects = ['A', 'B', 'C', 'D']
costs = {'A': 10, 'B': 15, 'C': 12, 'D': 8}
values = {'A': 20, 'B': 30, 'C': 25, 'D': 15}
budget = 30

prob = LpProblem("Project_Selection", LpMaximize)

# Binary variables: select[p] = 1 if project p is selected
select = LpVariable.dicts("select", projects, cat='Binary')

# Objective: maximize total value
prob += lpSum([values[p] * select[p] for p in projects])

# Budget constraint
prob += lpSum([costs[p] * select[p] for p in projects]) <= budget

# Solve
prob.solve()

# Display selected projects
print("Selected projects:")
for p in projects:
    if value(select[p]) == 1:
        print(f"  {p}: value={values[p]}, cost={costs[p]}")
print(f"Total value: {value(prob.objective)}")

Basic Example: Multi-Index Variables#

from pulp import *

# Transportation problem: ship products from factories to warehouses
factories = ['F1', 'F2', 'F3']
warehouses = ['W1', 'W2', 'W3', 'W4']

supply = {'F1': 100, 'F2': 150, 'F3': 120}
demand = {'W1': 80, 'W2': 70, 'W3': 90, 'W4': 60}
cost = {  # cost[factory, warehouse]
    ('F1', 'W1'): 2, ('F1', 'W2'): 4, ('F1', 'W3'): 5, ('F1', 'W4'): 3,
    ('F2', 'W1'): 3, ('F2', 'W2'): 1, ('F2', 'W3'): 3, ('F2', 'W4'): 2,
    ('F3', 'W1'): 5, ('F3', 'W2'): 4, ('F3', 'W3'): 2, ('F3', 'W4'): 3,
}

prob = LpProblem("Transportation", LpMinimize)

# Variables: flow[f, w] = amount shipped from factory f to warehouse w
routes = [(f, w) for f in factories for w in warehouses]
flow = LpVariable.dicts("flow", routes, lowBound=0)

# Objective: minimize total transportation cost
prob += lpSum([cost[f, w] * flow[f, w] for (f, w) in routes])

# Supply constraints: don't exceed factory capacity
for f in factories:
    prob += lpSum([flow[f, w] for w in warehouses]) <= supply[f]

# Demand constraints: meet warehouse demand
for w in warehouses:
    prob += lpSum([flow[f, w] for f in factories]) >= demand[w]

# Solve
prob.solve()

print(f"Total cost: {value(prob.objective)}")
for (f, w) in routes:
    if value(flow[f, w]) > 0:
        print(f"Ship {value(flow[f, w])} from {f} to {w}")

API Design Philosophy#

Algebraic modeling with Python syntax: Use Python operators (+, -, *, <=, >=, ==) naturally.

Advantages:

• Simple to learn: Minimal concepts (LpProblem, LpVariable, constraints)
• Pythonic: Uses Python operators, dictionaries, list comprehensions
• Bundled solver: CBC included, works immediately
• Lightweight: Focused API, not trying to do everything
• Readable: Models look like mathematical formulations

Disadvantages:

• Limited to LP/MILP: No nonlinear, no constraint programming
• Less flexible than Pyomo: Fewer solver options, less extensible
• Performance: For maximum MILP performance, commercial solvers (Gurobi) faster

Solver Backends#

PuLP interfaces with multiple LP/MILP solvers:

Bundled (Included):#

• CBC (COIN-OR Branch and Cut): Default, good quality open-source MILP solver

External (if installed):#

• GLPK: GNU Linear Programming Kit
• HiGHS: Modern high-performance LP/MIP solver
• Gurobi: Commercial (fast, requires license)
• CPLEX: IBM commercial (requires license)
• XPRESS: FICO commercial (requires license)
• SCIP: Open-source MILP/MINLP
• MOSEK: Commercial conic solver

Specify solver: prob.solve(GLPK()), prob.solve(GUROBI()), etc.

Install solvers via extras: pip install pulp[gurobi], pip install pulp[highs]

When to Use PuLP#

✅ Good fit when:#

• LP or MILP problems: Linear objective and constraints with integer variables
• Getting started: Learning optimization, simple problems
• Rapid prototyping: Quick experiments, proof of concept
• No external dependencies wanted: CBC bundled, works immediately
• Small to medium problems: Thousands of variables
• Occasional optimization: Not building optimization-heavy systems

❌ Not suitable when:#

• Nonlinear problems: No NLP support → use scipy, Pyomo, GEKKO
• Constraint programming: Scheduling, sequencing → use OR-Tools CP-SAT
• Convex optimization: SOCP, SDP → use CVXPY
• Need many solvers: Pyomo supports 20+ solvers, PuLP has fewer
• Very large MILP: For maximum performance, use Gurobi/CPLEX directly

Community and Maturity#

PuLP is a mature, stable library with active maintenance. Part of respected COIN-OR project.

Key Findings from Research#

1. Easiest entry point: Simplest API among algebraic modeling languages. Good for teaching and learning.

2. CBC bundled: Major usability advantage. No solver installation required. CBC is respectable open-source MILP solver.

3. Version 3.3.0 (Sept 2025): Recent release shows active maintenance. Requires Python 3.9+.

4. HiGHS support: Can use HiGHS (the rising star open-source solver adopted by SciPy and MATLAB) via pip install pulp[highs].

5. COIN-OR ecosystem: Part of mature operations research ecosystem (CBC, CLP, GLPK all COIN-OR projects).

Comparison with Alternatives#

Example Use Cases (Generic)#

Blending Problem#

Mix ingredients to meet requirements at minimum cost.

from pulp import *

# Blend ingredients to meet nutrient requirements
ingredients = ['A', 'B', 'C']
nutrients = ['protein', 'fiber', 'fat']

# Cost per unit
cost = {'A': 2.0, 'B': 3.5, 'C': 1.8}

# Nutrient content per unit
content = {
    ('A', 'protein'): 0.10, ('A', 'fiber'): 0.05, ('A', 'fat'): 0.02,
    ('B', 'protein'): 0.15, ('B', 'fiber'): 0.08, ('B', 'fat'): 0.03,
    ('C', 'protein'): 0.08, ('C', 'fiber'): 0.12, ('C', 'fat'): 0.01,
}

# Minimum nutrient requirements
requirements = {'protein': 10, 'fiber': 6, 'fat': 2}

prob = LpProblem("Blending", LpMinimize)

# Variables: amount of each ingredient
x = LpVariable.dicts("ingredient", ingredients, lowBound=0)

# Objective: minimize cost
prob += lpSum([cost[i] * x[i] for i in ingredients])

# Nutrient constraints
for n in nutrients:
    prob += lpSum([content[i, n] * x[i] for i in ingredients]) >= requirements[n]

prob.solve()

print(f"Minimum cost: {value(prob.objective):.2f}")
for i in ingredients:
    print(f"Ingredient {i}: {value(x[i]):.2f} units")

Cutting Stock Problem#

Cut raw materials to required lengths minimizing waste.

from pulp import *

# Cut stock of length 100 into required pieces
stock_length = 100
required_pieces = {30: 5, 40: 3, 50: 2}  # {length: quantity}

# Generate cutting patterns (simplified)
patterns = [
    {30: 3, 40: 0, 50: 0},  # Three 30s
    {30: 1, 40: 1, 50: 0},  # One 30, one 40
    {30: 0, 40: 2, 50: 0},  # Two 40s
    {30: 0, 40: 0, 50: 2},  # Two 50s
    {30: 1, 40: 0, 50: 1},  # One 30, one 50
]

prob = LpProblem("Cutting_Stock", LpMinimize)

# Variables: number of times each pattern is used
use_pattern = LpVariable.dicts("pattern", range(len(patterns)), lowBound=0, cat='Integer')

# Objective: minimize number of stocks used
prob += lpSum([use_pattern[p] for p in range(len(patterns))])

# Meet demand for each piece length
for length, quantity in required_pieces.items():
    prob += lpSum([patterns[p][length] * use_pattern[p]
                   for p in range(len(patterns))]) >= quantity

prob.solve()

print(f"Minimum stocks needed: {value(prob.objective)}")
for p in range(len(patterns)):
    if value(use_pattern[p]) > 0:
        print(f"Pattern {p}: use {value(use_pattern[p])} times - {patterns[p]}")

Facility Location#

Decide which facilities to open and how to serve customers.

from pulp import *

# Same as earlier example in or-tools.md
# (Facility location with fixed costs and transportation)
# ... (code similar to OR-Tools example but using PuLP syntax)

Advanced Features#

Multiple Objectives (Lexicographic)#

Solve for primary objective, then optimize secondary subject to primary being optimal.

# First solve for primary objective
prob.solve()
primary_optimal = value(prob.objective)

# Add constraint that primary stays optimal
prob += primary_objective == primary_optimal

# Change to secondary objective
prob.sense = LpMinimize
prob.setObjective(secondary_objective)

prob.solve()

Constraint Names and Access#

# Named constraints for later reference
prob += (x + y <= 10, "capacity_constraint")

# Access constraint
for name, constraint in prob.constraints.items():
    print(f"{name}: {constraint}")

Variable Bounds and Categories#

# Different variable types
x = LpVariable("x", lowBound=0)                    # x >= 0
y = LpVariable("y", lowBound=0, upBound=10)        # 0 <= y <= 10
z = LpVariable("z", cat='Integer')                  # Integer
b = LpVariable("b", cat='Binary')                   # Binary (0 or 1)

Writing Problem to File#

# Export to LP format for inspection or use with external tools
prob.writeLP("model.lp")

Performance Tips#

1. Use CBC for general MILP: Bundled, good performance
2. Use HiGHS for large LP: Install with pip install pulp[highs], faster than CBC for LP
3. Use Gurobi for large MILP: If performance critical and license available
4. Formulation matters: Tight formulations solve faster than loose ones
5. Use appropriate variable types: Don’t use Integer if Continuous works

Solver Selection Example#

from pulp import *

prob = LpProblem("Example", LpMaximize)
# ... define variables, objective, constraints ...

# Try different solvers
# prob.solve(PULP_CBC_CMD())        # Default bundled CBC
# prob.solve(HiGHS_CMD())           # HiGHS (if installed)
# prob.solve(GUROBI_CMD())          # Gurobi (if licensed)
# prob.solve(CPLEX_CMD())           # CPLEX (if licensed)
# prob.solve(GLPK_CMD())            # GLPK (if installed)

# With options
prob.solve(PULP_CBC_CMD(msg=1, timeLimit=60))  # Show output, 60s limit

References#

• Official documentation: https://coin-or.github.io/pulp/
• GitHub repository: https://github.com/coin-or/pulp
• COIN-OR: https://www.coin-or.org/
• Tutorials: https://coin-or.github.io/pulp/CaseStudies/index.html

Summary#

PuLP is the simplest path to LP/MILP optimization in Python:

• Easy to learn and use
• CBC solver bundled (works immediately)
• Clean algebraic modeling interface
• Mature and stable (COIN-OR project)
• Good for small to medium problems

Choose PuLP when:

• Getting started with optimization
• LP or MILP problems (no nonlinear)
• Want bundled solver (no installation hassles)
• Rapid prototyping
• Teaching or learning

Look elsewhere for:

• Nonlinear problems (scipy, Pyomo, GEKKO)
• Constraint programming (OR-Tools)
• Convex optimization with verification (CVXPY)
• Maximum solver flexibility (Pyomo’s 20+ backends)
• Large-scale performance (direct Gurobi/CPLEX API)

If you’re solving LP/MILP and want the easiest path, start with PuLP.

Pyomo: Python Optimization Modeling Objects#

Overview#

Pyomo is a Python-based open-source software package for formulating and analyzing optimization models. It is an algebraic modeling language (AML) that allows you to express optimization problems using mathematical notation, similar to AMPL or GAMS.

Key finding: With 2,127 GitHub stars and 123,641 weekly downloads, Pyomo is the leading Python algebraic modeling language, offering interfaces to 20+ solver backends.

Installation#

# Basic installation
pip install pyomo

# With optional solver interfaces
pip install pyomo[optional]

Note: Pyomo is a modeling language, not a solver. You must install solvers separately (GLPK, CBC, IPOPT, Gurobi, etc.).

Problem Types Supported#

Pyomo supports virtually all optimization problem types through its extensible architecture:

Pyomo’s strength is solver flexibility: same model code can be solved with different backends.

Basic Example: Linear Programming#

from pyomo.environ import *

# Create a concrete model
model = ConcreteModel()

# Decision variables
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)

# Objective: maximize x + 2y
model.obj = Objective(expr=model.x + 2*model.y, sense=maximize)

# Constraints
model.con1 = Constraint(expr=2*model.x + model.y <= 20)
model.con2 = Constraint(expr=model.x + 2*model.y <= 20)

# Solve with GLPK (open-source LP solver)
solver = SolverFactory('glpk')
result = solver.solve(model, tee=True)  # tee=True shows solver output

# Display results
print(f"x = {value(model.x)}")
print(f"y = {value(model.y)}")
print(f"Objective = {value(model.obj)}")

Basic Example: Mixed-Integer Programming#

from pyomo.environ import *

model = ConcreteModel()

# Integer decision variables
model.x = Var(domain=NonNegativeIntegers)
model.y = Var(domain=NonNegativeIntegers)

# Objective
model.obj = Objective(expr=model.x + 2*model.y, sense=maximize)

# Constraints
model.con1 = Constraint(expr=2*model.x + model.y <= 20)
model.con2 = Constraint(expr=model.x + 2*model.y <= 20)

# Solve with CBC (open-source MILP solver)
solver = SolverFactory('cbc')
result = solver.solve(model)

print(f"x = {value(model.x)}")  # Integer value
print(f"y = {value(model.y)}")  # Integer value
print(f"Objective = {value(model.obj)}")

Basic Example: Nonlinear Programming#

from pyomo.environ import *

model = ConcreteModel()

# Variables with bounds
model.x = Var(bounds=(0, 10), initialize=5)
model.y = Var(bounds=(0, 10), initialize=5)

# Nonlinear objective
model.obj = Objective(expr=model.x**2 + model.y**2)

# Nonlinear constraint
model.con1 = Constraint(expr=model.x**2 + model.y >= 4)

# Solve with IPOPT (open-source NLP solver)
solver = SolverFactory('ipopt')
result = solver.solve(model, tee=False)

print(f"x = {value(model.x)}")
print(f"y = {value(model.y)}")
print(f"Objective = {value(model.obj)}")

Basic Example: Indexed Variables and Constraints#

Pyomo excels at modeling problems with sets and indices:

from pyomo.environ import *

model = ConcreteModel()

# Sets
model.PROJECTS = Set(initialize=['A', 'B', 'C', 'D'])

# Parameters (data)
cost = {'A': 10, 'B': 15, 'C': 12, 'D': 8}
return_val = {'A': 20, 'B': 30, 'C': 25, 'D': 15}
model.budget = Param(initialize=30)

# Indexed binary variables: select[p] = 1 if project p is selected
model.select = Var(model.PROJECTS, domain=Binary)

# Objective: maximize total return
model.obj = Objective(
    expr=sum(return_val[p] * model.select[p] for p in model.PROJECTS),
    sense=maximize
)

# Budget constraint
model.budget_con = Constraint(
    expr=sum(cost[p] * model.select[p] for p in model.PROJECTS) <= model.budget
)

# Solve
solver = SolverFactory('cbc')
solver.solve(model)

# Display selected projects
print("Selected projects:")
for p in model.PROJECTS:
    if value(model.select[p]) > 0.5:  # Binary variable
        print(f"  {p}: return = {return_val[p]}, cost = {cost[p]}")

API Design Philosophy#

Algebraic modeling language: Express models using mathematical notation close to textbook formulations.

Advantages:

• Readable models (looks like mathematics)
• Separation of model and data
• Reusable model structure
• Easy to modify and extend
• Solver-agnostic (change solver with one line)
• Structured indexing (sets, parameters, indexed variables)

Disadvantages:

• Learning curve (new syntax and concepts)
• More verbose than direct solver APIs for simple problems
• Overhead for very small problems
• Requires understanding of both Pyomo syntax and optimization theory

Solver Backends (20+ supported)#

Pyomo interfaces with external solvers via file-based or direct API communication.

Open-Source Linear/Integer#

• GLPK: GNU Linear Programming Kit
• CBC: COIN-OR Branch and Cut
• HiGHS: High-performance LP/MIP solver

Open-Source Nonlinear#

• IPOPT: Interior Point OPTimizer (large-scale NLP)
• Couenne: Convex Over and Under ENvelopes for NLP (MINLP)
• BONMIN: Basic Open-source Nonlinear Mixed INteger (MINLP)

Commercial#

• Gurobi: Leading commercial MILP/QP solver
• CPLEX: IBM commercial solver
• MOSEK: Conic optimization specialist
• XPRESS: FICO solver
• KNITRO: Nonlinear optimization
• BARON: Global MINLP solver

Specialized#

• Mindtpy: Pyomo’s built-in MINLP solver
• SCIP: Open-source MILP/MINLP
• GDPopt: Generalized disjunctive programming

Use SolverFactory('solver_name') to select solver.

When to Use Pyomo#

✅ Good fit when:#

• Complex models: Many constraints, indexed variables, reusable structure
• Solver flexibility: Want to try multiple solvers (CBC, Gurobi, IPOPT) on same model
• Academic research: Need to document model structure, try different formulations
• Data-driven models: Separate model logic from data, load data from files/databases
• Mixed problem types: LP, MILP, NLP in same project
• Advanced features: Stochastic programming, disjunctive programming, DAE systems
• Model archival: Mathematical formulation readable years later

❌ Not suitable when:#

• Simple one-off problems: scipy.optimize or PuLP simpler
• Constraint programming: OR-Tools CP-SAT better for scheduling
• Convex optimization: CVXPY provides DCP verification
• Maximum performance: Direct solver APIs (gurobipy) slightly faster
• No solver installation: Pyomo requires external solvers (unlike OR-Tools)

Community and Maturity#

Pyomo is the academic standard for optimization modeling in Python. Used extensively in research and teaching.

Key Findings from Research#

1. Ecosystem leader: 20+ solver interfaces make Pyomo the most flexible Python modeling language for solver experimentation.

2. Active development: Regular releases (v6.9.4 in 2025), new features like logic-based discrete steepest descent in GDPopt.

3. Extensions: Pyomo.DAE (differential-algebraic equations), PySP (stochastic programming), GDP (disjunctive programming) enable advanced problem types.

4. Educational adoption: Many universities use Pyomo for teaching optimization. Textbook “Pyomo — Optimization Modeling in Python” available.

5. Academic roots: Developed at Sandia National Laboratories, widely used in energy systems, chemical engineering, operations research.

Comparison with Alternatives#

Example Use Cases (Generic)#

Production Planning with Time Periods#

Determine production quantities over time to minimize cost while meeting demand.

from pyomo.environ import *

model = ConcreteModel()

# Sets
model.PRODUCTS = Set(initialize=['A', 'B', 'C'])
model.PERIODS = RangeSet(1, 12)  # 12 time periods

# Parameters (data)
model.demand = Param(model.PRODUCTS, model.PERIODS, initialize=demand_data)
model.production_cost = Param(model.PRODUCTS, initialize=prod_cost_data)
model.inventory_cost = Param(model.PRODUCTS, initialize=inv_cost_data)
model.capacity = Param(model.PERIODS, initialize=capacity_data)

# Variables
model.produce = Var(model.PRODUCTS, model.PERIODS, domain=NonNegativeReals)
model.inventory = Var(model.PRODUCTS, model.PERIODS, domain=NonNegativeReals)

# Objective: minimize total cost
def total_cost_rule(m):
    return sum(m.production_cost[p] * m.produce[p, t] +
               m.inventory_cost[p] * m.inventory[p, t]
               for p in m.PRODUCTS for t in m.PERIODS)
model.obj = Objective(rule=total_cost_rule, sense=minimize)

# Inventory balance constraints
def inventory_balance_rule(m, p, t):
    if t == 1:
        return m.inventory[p, t] == m.produce[p, t] - m.demand[p, t]
    else:
        return m.inventory[p, t] == m.inventory[p, t-1] + m.produce[p, t] - m.demand[p, t]
model.inventory_balance = Constraint(model.PRODUCTS, model.PERIODS, rule=inventory_balance_rule)

# Capacity constraint
def capacity_rule(m, t):
    return sum(m.produce[p, t] for p in m.PRODUCTS) <= m.capacity[t]
model.capacity_con = Constraint(model.PERIODS, rule=capacity_rule)

solver = SolverFactory('glpk')
solver.solve(model)

Network Design with Fixed Costs#

Select facilities to open (binary decisions) and route flows to minimize total cost.

from pyomo.environ import *

model = ConcreteModel()

model.FACILITIES = Set(initialize=['F1', 'F2', 'F3'])
model.CUSTOMERS = Set(initialize=['C1', 'C2', 'C3', 'C4'])

# Binary: open facility or not
model.open = Var(model.FACILITIES, domain=Binary)

# Continuous: flow from facility to customer
model.flow = Var(model.FACILITIES, model.CUSTOMERS, domain=NonNegativeReals)

# Objective: minimize fixed costs + variable costs
def total_cost_rule(m):
    fixed = sum(fixed_cost[f] * m.open[f] for f in m.FACILITIES)
    variable = sum(transport_cost[f, c] * m.flow[f, c]
                   for f in m.FACILITIES for c in m.CUSTOMERS)
    return fixed + variable
model.obj = Objective(rule=total_cost_rule, sense=minimize)

# Meet customer demand
def demand_rule(m, c):
    return sum(m.flow[f, c] for f in m.FACILITIES) >= demand[c]
model.demand_con = Constraint(model.CUSTOMERS, rule=demand_rule)

# Facility capacity (only if open)
def capacity_rule(m, f):
    return sum(m.flow[f, c] for c in m.CUSTOMERS) <= capacity[f] * m.open[f]
model.capacity_con = Constraint(model.FACILITIES, rule=capacity_rule)

solver = SolverFactory('cbc')
solver.solve(model)

Performance Characteristics#

Pyomo adds modest overhead for:

• Model construction (building expression trees)
• Communication with solver (file-based or API)

For large models (10,000+ variables), overhead typically <10% of total solve time. Solver performance dominates.

Tips for performance:

• Use indexed variables/constraints instead of loops
• Choose fast solvers (Gurobi > CBC for MILP)
• Consider direct API interfaces (pyomo-gurobi-direct) for large models

Advanced Features#

Abstract Models#

Separate model structure from data:

model = AbstractModel()
# Define structure with undefined parameters
model.load('data.dat')  # Load data later

Persistent Solvers#

Reuse solver instance across multiple solves (incremental changes):

solver = SolverFactory('gurobi', solver_io='python')
solver.set_instance(model)
solver.solve()
# Modify model
solver.solve()  # Incremental solve

Suffixes#

Access solver-specific information (dual values, reduced costs):

model.dual = Suffix(direction=Suffix.IMPORT)
solver.solve(model)
print(model.dual[model.con1])  # Dual value of constraint

References#

• Official documentation: https://pyomo.readthedocs.io
• GitHub: https://github.com/Pyomo/pyomo (2.1k stars)
• Textbook: Bynum et al., “Pyomo — Optimization Modeling in Python”, Springer
• Website: http://www.pyomo.org

Summary#

Pyomo is the most comprehensive Python optimization modeling language:

• 20+ solver backends (most flexible)
• All problem types (LP, MILP, NLP, MINLP, stochastic, dynamic)
• Algebraic modeling for complex, structured problems
• Academic standard for optimization research

Choose Pyomo when:

• You need solver flexibility (try Gurobi, CBC, IPOPT on same model)
• Model complexity justifies algebraic modeling overhead
• Academic research or teaching
• Advanced features (stochastic, disjunctive, DAE)

Look elsewhere for:

• Simple one-off problems (scipy.optimize, PuLP)
• Convex with verification (CVXPY)
• Constraint programming (OR-Tools)
• Maximum performance (direct APIs)

S1-Rapid: Library Selection Recommendations#

Executive Summary#

Based on research of Python optimization libraries (PyPI statistics, GitHub activity, documentation quality, and community feedback), selection depends primarily on problem type, followed by solver flexibility needs and simplicity requirements.

Key insight: Problem type (LP, MILP, NLP, convex) is the dominant selection criterion, not application domain.

Quick Decision Tree#

What problem type?
│
├─ Linear/Integer Programming (LP/MILP)?
│  ├─ Simple, one-off problem? → PuLP (easiest)
│  ├─ Need solver flexibility? → Pyomo (20+ solvers)
│  ├─ Production, constraint programming? → OR-Tools (CP-SAT)
│  └─ Already using scipy? → scipy.optimize.linprog (HiGHS backend)
│
├─ Convex Optimization?
│  └─ CVXPY (DCP verification, specialized solvers)
│
├─ Nonlinear Programming (NLP)?
│  ├─ Small-medium scale? → scipy.optimize (BFGS, SLSQP)
│  ├─ Large scale or need solvers? → Pyomo + IPOPT
│  └─ Dynamic optimization / DAE? → GEKKO
│
├─ Multi-Objective?
│  └─ pymoo (Pareto frontier, genetic algorithms)
│
└─ Constraint Programming / Scheduling?
   └─ OR-Tools CP-SAT (best-in-class)

Detailed Recommendations by Scenario#

Scenario 1: Getting Started / Learning#

Recommendation: PuLP (if LP/MILP) or scipy.optimize (if continuous)

Rationale:

• PuLP: Simplest algebraic modeling for LP/MILP, CBC bundled
• scipy.optimize: Zero installation (part of scipy), good for continuous problems
• Both have gentle learning curves and excellent documentation

Example: Student learning optimization, analyst solving occasional problems

Scenario 2: Rapid Prototyping#

Recommendation: scipy.optimize (continuous) or PuLP (discrete)

Rationale:

• Fast to code, minimal boilerplate
• Good for experiments and proof-of-concept
• Can graduate to more powerful tools later

Example: Researcher testing optimization approach, data scientist exploring models

Scenario 3: Production LP/MILP#

Recommendation: OR-Tools or Pyomo + HiGHS/Gurobi

Rationale:

• OR-Tools: Google-maintained, bundled solvers, excellent performance (CP-SAT), production-ready
• Pyomo + HiGHS: Open-source, competitive performance, flexible
• Pyomo + Gurobi: Maximum performance if license available

Example: Logistics optimization service, resource allocation system, production planning application

Scenario 4: Convex Optimization#

Recommendation: CVXPY

Rationale:

• Automatic convexity verification (DCP analysis)
• Specialized solvers for convex problems
• Standard tool for machine learning, portfolio optimization, signal processing

Example: Portfolio optimization, Lasso/ridge regression, SVM training, control systems

Scenario 5: Nonlinear Programming#

Recommendation: scipy.optimize (small-medium) or Pyomo + IPOPT (large)

Rationale:

• scipy.optimize: Built-in, good algorithms (BFGS, trust-region), handles constraints
• Pyomo + IPOPT: Large-scale NLP, more solver options, algebraic modeling

Example: Parameter estimation, engineering design optimization, calibration against simulation

Scenario 6: Scheduling / Combinatorial#

Recommendation: OR-Tools CP-SAT

Rationale:

• Best-in-class constraint programming solver
• Specialized for scheduling, sequencing, assignment
• MiniZinc Challenge winner

Example: Job shop scheduling, employee rostering, task assignment, timetabling

Scenario 7: Vehicle Routing#

Recommendation: OR-Tools routing module

Rationale:

• Specialized routing solvers
• Rich constraint modeling (time windows, capacities, precedence)
• Tuned metaheuristics

Example: Delivery routing, technician dispatch, field service optimization

Scenario 8: Academic Research / Multi-Solver Comparison#

Recommendation: Pyomo

Rationale:

• 20+ solver backends (most flexible)
• Same model code across different solvers
• Extensive features (stochastic, disjunctive, DAE)
• Academic standard

Example: Comparing MILP solvers, testing problem formulations, operations research research

Scenario 9: Maximum MILP Performance#

Recommendation: Gurobi (direct API) or CPLEX (direct API)

Rationale:

• Fastest commercial solvers for large-scale MILP
• Direct API (gurobipy, docplex) avoids modeling language overhead
• Worth cost for time-critical, large-scale problems

Example: Large-scale logistics (100k+ variables), real-time optimization with tight deadlines

Scenario 10: Multi-Objective Optimization#

Recommendation: pymoo

Rationale:

• Only library specialized for multi-objective
• Generates Pareto frontier
• Modern evolutionary algorithms (NSGA-II, NSGA-III, MOEAD)

Example: Design optimization with competing objectives, trade-off analysis

Scenario 11: Dynamic Optimization / Optimal Control#

Recommendation: GEKKO or Pyomo.DAE

Rationale:

• GEKKO: Specialized for dynamic optimization, MPC, DAE systems
• Pyomo.DAE: Integrates with Pyomo ecosystem

Example: Process control, trajectory optimization, model predictive control

Scenario 12: Already Using SciPy/NumPy#

Recommendation: scipy.optimize (first), consider CVXPY if convex

Rationale:

• Zero additional dependencies
• Integrates with existing scipy/numpy code
• Good default for scientific computing workflows

Example: Scientific computing pipelines, data analysis notebooks

Library Comparison Matrix#

Problem Type → Library Mapping#

Solver Selection Guidance#

Open-Source Solvers (Free)#

Linear Programming:

• HiGHS: Rising star, adopted by SciPy and MATLAB. Fast, MIT license.
• GLPK: Mature, widely available. Slower than HiGHS.
• CLP: COIN-OR, good quality.

Mixed-Integer:

• SCIP: Fastest academic solver. Now Apache 2.0 (was academic-only).
• CBC: COIN-OR, bundled with PuLP. Solid performance.
• HiGHS: Also does MILP (not just LP).

Nonlinear:

• IPOPT: Interior point, large-scale. EPL license.
• Bonmin, Couenne: MINLP solvers. Academic use.

Constraint Programming:

• OR-Tools CP-SAT: Best-in-class. Apache 2.0.

Convex:

• Clarabel: Modern, Rust-based. Default in CVXPY.
• SCS: Robust, medium accuracy.
• OSQP: Fast for QP.

Commercial Solvers (Licensed)#

When to consider:

• Large-scale MILP (100k+ variables)
• Time-critical applications
• Maximum performance needed
• Budget available

Top choices:

• Gurobi: Fastest overall, excellent support, free for academics
• CPLEX: IBM, strong on high-dimensionality, free for academics
• MOSEK: Best for conic (SOCP, SDP), free for academics

Cost vs benefit: For small-medium problems, open-source competitive. For very large or time-critical, commercial worth cost.

Common Migration Paths#

1. Learning → Production#

Path: PuLP → OR-Tools or Pyomo + HiGHS

Rationale: Start simple with PuLP, graduate to production-grade tools as needs grow.

2. Prototype → Scale#

Path: scipy.optimize → Pyomo + IPOPT → Pyomo + commercial NLP solver

Rationale: Prototype with scipy, add algebraic modeling for complexity, scale to commercial for performance.

3. Open-Source → Commercial#

Path: Pyomo + CBC → Pyomo + Gurobi (same model code!)

Rationale: Develop with open-source, deploy with commercial for performance. Pyomo makes this seamless.

4. General → Specialized#

Path: Pyomo → CVXPY (if convex) or OR-Tools CP-SAT (if scheduling)

Rationale: Recognize problem structure, use specialized tool.

Anti-Recommendations#

❌ Don’t use scipy.optimize for:#

• Integer variables (no MILP support)
• Large-scale LP (use HiGHS via Pyomo or scipy.optimize.linprog)

❌ Don’t use PuLP for:#

• Nonlinear problems (no NLP support)
• Constraint programming (no CP support)

❌ Don’t use CVXPY for:#

• Non-convex problems (DCP will reject)
• Integer programming (limited support)

❌ Don’t use OR-Tools for:#

• Nonlinear programming (no NLP support)
• If algebraic modeling strongly preferred (Pyomo better)

❌ Don’t use Pyomo for:#

• Simple one-off problems (overhead not worth it)
• Constraint programming (OR-Tools CP-SAT better)

❌ Don’t use commercial solvers for:#

• Small problems (overkill)
• Budget-constrained projects (open-source competitive)
• Academic use without free license

Final Recommendations Summary#

Default Starting Points#

1. LP/MILP: Start with PuLP (easiest)
2. Convex: Use CVXPY (DCP verification)
3. Continuous NLP: Use scipy.optimize (built-in)
4. Scheduling: Use OR-Tools CP-SAT (best-in-class)
5. Multi-objective: Use pymoo (specialized)

When to Graduate to More Powerful Tools#

• PuLP → Pyomo when need solver flexibility or complex models
• scipy.optimize → Pyomo + IPOPT when NLP gets large
• Open-source → Commercial when performance critical (Gurobi, CPLEX)

Production Deployments#

• OR-Tools: Safest choice (Google-backed, bundled solvers, proven)
• Pyomo + HiGHS/Gurobi: Strong alternative with solver flexibility
• CVXPY: If convex problems (finance, ML, control)

Academic Research#

• Pyomo: Standard for operations research (20+ solvers, extensive features)
• CVXPY: Standard for convex optimization (Stanford, widely cited)
• OR-Tools: For constraint programming research

Surprising Findings from Research#

1. HiGHS adoption: SciPy and MATLAB both chose HiGHS as default LP/MIP solver. Strong signal of quality.

2. SCIP now Apache 2.0: Previously academic-only, now fully open source. Major development for open-source MILP.

3. Gurobi/CPLEX withdrew from Mittelmann benchmarks: Reduces transparency. Open-source performance harder to compare.

4. OR-Tools CP-SAT dominance: Google’s CP-SAT wins MiniZinc Challenge medals, competitive with commercial CP solvers.

5. CVXPY’s DCP is unique: Only Python library with automatic convexity verification. Killer feature for convex problems.

6. Pyomo’s 20+ solvers: Far exceeds other modeling languages in flexibility.

7. PuLP’s simplicity: Still best entry point after 15+ years.

References#

This recommendation is based on:

• Library download statistics (PyPI, conda)
• GitHub metrics (stars, activity, contributors)
• Documentation quality assessment
• Academic publications (JMLR, IEEE Access, Nature Methods)
• Mittelmann benchmark analysis
• Community feedback (Stack Overflow, forums)
• Solver adoption patterns (SciPy, MATLAB)

Conclusion#

Problem type determines library choice. Once you know your problem type:

• LP/MILP: PuLP (simple) or OR-Tools/Pyomo (production)
• Convex: CVXPY
• NLP: scipy.optimize (small) or Pyomo (large)
• CP: OR-Tools CP-SAT
• Multi-objective: pymoo

Start with the simplest tool that fits your problem. Graduate to more powerful tools as needs grow. Pyomo provides migration path from open-source to commercial solvers with same model code.

The Python optimization ecosystem in 2025 is mature, with excellent open-source options competitive with commercial tools for many use cases.

scipy.optimize: Python’s Built-In Optimization#

Overview#

scipy.optimize is part of SciPy, the foundational scientific computing library for Python. It provides optimization algorithms for continuous problems, both constrained and unconstrained. Part of the scipy package which has millions of downloads per year (>13M from PyPI in 2019, likely much higher now as part of scientific Python stack).

Key insight: scipy.optimize is ubiquitous but limited to smaller-scale problems and lacks mixed-integer programming capability.

Installation#

pip install scipy

scipy is typically pre-installed in scientific Python distributions (Anaconda, etc.) and has no external dependencies beyond NumPy.

Problem Types Supported#

Basic Example: Unconstrained Optimization#

import numpy as np
from scipy.optimize import minimize

# Define Rosenbrock function (classic test problem)
def rosenbrock(x):
    return sum(100.0 * (x[1:] - x[:-1]**2)**2 + (1 - x[:-1])**2)

# Initial guess
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])

# Optimize with BFGS (gradient-based)
result = minimize(rosenbrock, x0, method='BFGS')

print(f"Optimal solution: {result.x}")
print(f"Optimal value: {result.fun}")
print(f"Success: {result.success}")

Basic Example: Constrained Optimization#

from scipy.optimize import minimize, LinearConstraint, NonlinearConstraint

# Minimize x^2 + y^2 subject to x + y >= 1, x >= 0, y >= 0
def objective(x):
    return x[0]**2 + x[1]**2

# Linear constraint: x + y >= 1  →  -x - y <= -1
linear_constraint = LinearConstraint([[1, 1]], [1], [np.inf])

# Bounds: x >= 0, y >= 0
bounds = [(0, None), (0, None)]

# Optimize with SLSQP (Sequential Least Squares Programming)
result = minimize(objective, x0=[0, 0], method='SLSQP',
                  constraints=linear_constraint, bounds=bounds)

print(f"Optimal solution: {result.x}")  # Should be approximately [0.5, 0.5]
print(f"Optimal value: {result.fun}")   # Should be approximately 0.5

Basic Example: Linear Programming#

from scipy.optimize import linprog

# Minimize: c^T x
c = [-1, -2]  # Coefficients (negative because linprog minimizes)

# Subject to: A_ub @ x <= b_ub
A_ub = [[2, 1],   # 2x + y <= 20
        [1, 2]]   # x + 2y <= 20

b_ub = [20, 20]

# Bounds: x >= 0, y >= 0
bounds = [(0, None), (0, None)]

# Solve with HiGHS (default since SciPy 1.6.0)
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

print(f"Optimal solution: {result.x}")
print(f"Optimal value: {-result.fun}")  # Negate because we minimized negative

Available Optimization Methods#

Unconstrained Minimization (minimize)#

Constrained Minimization#

Global Optimization#

Linear Programming#

API Design Philosophy#

Functional interface: Pass objective function and constraints as Python functions. Solver calls them during optimization.

Advantages:

• Quick to get started
• Flexible (any Python code in objective)
• No need to learn modeling language

Disadvantages:

• No symbolic differentiation (relies on numerical gradients unless you provide analytical)
• Difficult to inspect or modify model structure
• Less efficient for large-scale problems
• Cannot export model to other solvers

Solver Backends#

scipy.optimize is both library and solver:

• Implements algorithms directly in Python/Cython (BFGS, Nelder-Mead, etc.)
• Wraps external solvers (HiGHS for LP since 1.6.0)

Cannot interface with external MILP or NLP solvers (Gurobi, CPLEX, IPOPT). For multi-solver support, use Pyomo or PuLP.

When to Use scipy.optimize#

✅ Good fit when:#

• Already using scipy/numpy: Zero additional dependencies
• Small to medium scale: Hundreds to thousands of variables
• Continuous variables: No integer decisions
• Rapid prototyping: Quick experiments, research code
• Function-based objective: Simulation output, machine learning model, complex calculation
• Scientific computing context: Parameter estimation, curve fitting, model calibration

❌ Not suitable when:#

• Integer variables required: No MILP capability → use PuLP, Pyomo, OR-Tools
• Large-scale LP/MILP: Thousands of constraints → use dedicated solvers
• Need solver flexibility: Cannot swap solvers (Gurobi vs CBC) → use Pyomo
• Complex model structure: Many constraints, reusable model → use algebraic modeling
• Production performance critical: Direct solver APIs (gurobipy) may be faster

Community and Maturity#

scipy is the foundation of scientific Python. Extremely stable, well-tested, and maintained.

Key Findings from Research#

1. HiGHS adoption (2021): SciPy 1.6.0 switched to HiGHS for linprog, providing competitive LP performance with open-source solver.

2. No MILP support: This is the biggest limitation. For integer variables, must use other libraries.

3. Gradient computation: If you don’t provide analytical gradients, scipy uses finite differences (slow, inaccurate). Consider automatic differentiation (JAX, autograd).

4. Trust-region methods: trust-constr (added 2018) is modern, robust method for constrained optimization. Recommended over older SLSQP for difficult problems.

5. Global optimization limited: Metaheuristics (differential_evolution) can escape local optima but no guarantee of global optimum. For convex problems, local methods sufficient.

Comparison with Alternatives#

Example Use Cases (Generic)#

Parameter Estimation#

Fit model parameters to data by minimizing prediction error.

def model_prediction(params, inputs):
    # Your model here
    return predictions

def objective(params):
    predictions = model_prediction(params, training_data)
    return np.sum((predictions - observations)**2)

result = minimize(objective, initial_params, method='L-BFGS-B')

Optimal Resource Allocation (Continuous)#

Allocate continuous quantities (e.g., budget, energy) across options to maximize utility.

def utility(allocation):
    # Compute total utility from allocation
    return -np.sum(utility_functions(allocation))  # Negative for minimization

constraints = {'type': 'eq', 'fun': lambda x: np.sum(x) - total_budget}
bounds = [(0, None)] * n_options

result = minimize(utility, initial_allocation, method='SLSQP',
                  constraints=constraints, bounds=bounds)

Calibration Against Simulation#

Find input parameters that make simulation output match observed data.

def simulation_error(params):
    # Run simulation with params
    sim_output = run_simulation(params)
    return np.sum((sim_output - target_output)**2)

result = minimize(simulation_error, initial_params, method='Nelder-Mead')

References#

• SciPy documentation: https://docs.scipy.org/doc/scipy/reference/optimize.html
• SciPy 1.0 paper: Virtanen et al., Nature Methods 2020
• HiGHS integration: SciPy 1.6.0 release notes (2021)

Summary#

scipy.optimize is the default starting point for continuous optimization in Python. Excellent for:

• Research and prototyping
• Parameter estimation and curve fitting
• Small to medium problems
• Users already in scipy/numpy ecosystem

For production MILP, multi-solver support, or large-scale problems, consider Pyomo, PuLP, or OR-Tools instead.

S2-Comprehensive: Commercial vs Open-Source Solvers#

Performance Gap Analysis#

Historical Context (Pre-2024)#

Commercial solvers (Gurobi, CPLEX) were 5-10x faster than open-source on large MILP instances

Current Status (2025)#

Major shift: Gurobi and MindOpt withdrew from Mittelmann benchmarks (Aug/Dec 2024), reducing transparency

Open-Source Solvers#

HiGHS (MIT License)#

Adoption signal: Chosen as default LP/MIP solver by:

• SciPy 1.6.0+ (2021)
• MATLAB Optimization Toolbox

Performance: Competitive with commercial on LP, good on MILP Conclusion: Rising star, closing performance gap

SCIP (Apache 2.0, since v9.0)#

Status: Fastest academic MILP/MINLP solver Previous barrier: Academic-only license Now: Fully open-source (Apache 2.0) Performance: Competitive with commercial on many problems

IPOPT (EPL)#

Focus: Large-scale nonlinear programming Performance: Industry standard open-source NLP solver Limitation: Local optima only (non-convex problems)

CBC (EPL)#

Maturity: 15+ years Bundled: PuLP default Performance: 2-5x slower than SCIP, adequate for small-medium MILP

Commercial Solvers#

Gurobi#

Pricing: ~$10k-100k+/year (enterprise) Academic: Free Performance: Historically fastest MILP solver Benchmark status: Withdrew August 2024

CPLEX (IBM)#

Pricing: Similar to Gurobi Academic: Free Strengths: High-dimensionality problems, non-convex MIQP Benchmark status: Still participates

MOSEK#

Focus: Conic optimization (SOCP, SDP) Pricing: ~$5k-50k/year Academic: Free Assessment: Best commercial option for conic problems

KNITRO (Artelys)#

Focus: Nonlinear programming Performance: 2-10x faster than IPOPT on many NLP Pricing: $5k-50k/year

When Commercial Worth Cost#

Clear Justifications:#

1. Very large MILP: 100k+ variables, 10k+ integer variables
2. Time-critical production: Real-time optimization, tight SLAs
3. ROI analysis: If optimization saves $100k+/year, $10k solver cost justified
4. Professional support: Mission-critical applications need vendor support

Questionable:#

1. Small-medium problems: Open-source competitive
2. Research/exploration: Need flexibility more than maximum speed
3. Budget constraints: $10k-100k is significant

Academic vs Industry Access#

Academic Advantages:#

• Free commercial licenses: Gurobi, CPLEX, MOSEK free for academic use
• Full functionality: No feature limitations
• Ideal for research: Test multiple solvers

Industry Reality:#

• Must pay: Commercial licenses expensive
• Open-source competitive: HiGHS, SCIP often sufficient
• Hybrid approach: Develop with open-source, deploy with commercial if needed

Cost Comparison#

Performance Tiers (Based on Historical Data)#

LP#

1. Tier 1: HiGHS, Gurobi, CPLEX (gap narrowed)
2. Tier 2: GLOP, CLP
3. Tier 3: GLPK

MILP#

1. Tier 1: Gurobi, CPLEX (historically)
2. Tier 2: SCIP (best open-source)
3. Tier 3: HiGHS, CBC
4. Tier 4: GLPK

NLP#

1. Tier 1: KNITRO, SNOPT
2. Tier 2: IPOPT (excellent for open-source)

Convex Conic#

1. Tier 1: MOSEK
2. Tier 2: Clarabel, SCS, ECOS

Open-Source Advantages#

1. Cost: Free (obviously)
2. Transparency: Source code available
3. Community: Active development, bug reporting
4. Flexibility: Modify if needed
5. No licensing: No license servers, floating licenses, compliance issues

Commercial Advantages#

1. Performance: Still faster on largest instances
2. Support: Professional support contracts
3. Features: Advanced tuning, cloud, distributed solving
4. Testing: Extensive QA, stability guarantees

Trend: Closing Gap#

Evidence:

• HiGHS adoption by SciPy, MATLAB
• SCIP now Apache 2.0 (removed licensing barrier)
• OR-Tools CP-SAT award-winning (beats commercial CP)

Implication: Open-source increasingly viable for production

Recommendations#

Start with Open-Source:#

• HiGHS (LP/MILP)
• SCIP (MILP/MINLP)
• IPOPT (NLP)
• OR-Tools (CP, routing)
• CVXPY with open solvers (convex)

Upgrade to Commercial When:#

• Profiling shows solver is bottleneck
• Problem size requires maximum performance
• Budget justifies cost (ROI positive)
• Professional support needed

Migration Path:#

Use Pyomo for seamless open-source → commercial migration:

# Development
solver = SolverFactory('scip')

# Production (same model code!)
solver = SolverFactory('gurobi')

Licensing Considerations#

Open-Source Licenses:#

• MIT (HiGHS): Most permissive
• Apache 2.0 (SCIP, OR-Tools): Permissive, patent protection
• EPL (IPOPT, CBC): Weak copyleft
• GPL: Strong copyleft (few optimization solvers)

Commercial Licenses:#

• Named-user: Tied to individual
• Floating: Shared pool of licenses
• Cloud: Usage-based pricing
• Academic: Free but restricted to academic use

Key Findings#

1. Open-source competitive for many use cases: No longer automatically commercial
2. HiGHS emergence: Major development (SciPy, MATLAB adoption)
3. SCIP Apache 2.0: Removes licensing barrier for best open-source MILP
4. Benchmark withdrawal: Reduced transparency for commercial performance
5. Commercial still faster on huge MILP: But gap narrowed
6. CP-SAT excellence: OR-Tools competitive with commercial CP solvers
7. Convex exception: MOSEK significantly faster than open-source for conic

Decision Framework#

Problem size and performance requirements?
├─ Small-medium → Open-source sufficient
│  └─ HiGHS, SCIP, IPOPT, OR-Tools
│
├─ Large, but not time-critical → Try open-source first
│  └─ Profile before deciding on commercial
│
└─ Very large AND time-critical → Commercial justified
   └─ Gurobi, CPLEX (MILP), MOSEK (conic), Knitro (NLP)

Budget available?
├─ No → Open-source (HiGHS, SCIP excellent)
├─ Academic → Free commercial licenses available
└─ Industry → Evaluate ROI (does speed improvement justify cost?)

Conclusion#

2025 landscape: Open-source solvers increasingly viable for production. Start with open-source; upgrade to commercial only when profiling demonstrates clear need. The performance gap has narrowed significantly, especially for LP and medium-scale MILP.

Biggest winners: HiGHS (LP/MILP), SCIP (MILP/MINLP), OR-Tools (CP), IPOPT (NLP)

Commercial still justified for: Very large MILP, professional support needs, conic optimization (MOSEK)

S2-Comprehensive: Modeling Languages vs Direct APIs#

Overview#

Optimization problems can be expressed through modeling languages (algebraic, declarative) or direct solver APIs (programmatic). Each approach has trade-offs.

Modeling Language Paradigm#

Concept#

Express optimization problems using mathematical notation similar to textbook formulations.

Examples#

• Pyomo: Python algebraic modeling
• CVXPY: Convex optimization modeling
• PuLP: Simple LP/MILP modeling
• AMPL, GAMS: Commercial modeling languages
• JuMP: Julia modeling language

Characteristics#

Algebraic syntax:

# Pyomo example
model = ConcreteModel()
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)
model.obj = Objective(expr=model.x + 2*model.y, sense=maximize)
model.con1 = Constraint(expr=2*model.x + model.y <= 20)

Advantages:

1. Readable: Looks like mathematics
2. Solver-agnostic: Change solver with one line
3. Separation of model and data: Reusable structure
4. Structured indexing: Sets, parameters, indexed variables
5. Model introspection: Examine structure programmatically
6. Educational: Clear formulation for teaching

Disadvantages:

1. Learning curve: New syntax and concepts
2. Performance overhead: 5-10% (typically)
3. Abstraction layer: Less control over solver internals
4. Installation: Modeling language + solver

When to Use#

• Complex models with many constraints
• Need solver flexibility
• Model reuse and maintenance
• Academic research
• Teaching optimization

Direct Solver API Paradigm#

Concept#

Call solver functions directly to build and solve models.

Examples#

• gurobipy: Gurobi Python API
• docplex: CPLEX Python API
• scipy.optimize: Function-based optimization
• OR-Tools: Object-oriented APIs

Characteristics#

Programmatic syntax:

# Gurobi direct API example
import gurobipy as gp

model = gp.Model()
x = model.addVar(name="x")
y = model.addVar(name="y")
model.setObjective(x + 2*y, gp.GRB.MAXIMIZE)
model.addConstr(2*x + y <= 20)
model.optimize()

Advantages:

1. Performance: Minimal overhead
2. Fine control: Access solver-specific features
3. Simpler installation: Just the solver
4. Incremental building: Add/remove constraints efficiently

Disadvantages:

1. Solver lock-in: Code tied to specific solver
2. Less readable: More verbose than algebraic
3. No solver abstraction: Can’t easily swap solvers

When to Use#

• Maximum performance needed
• Solver-specific features required
• Incremental model building (add/remove constraints frequently)
• Simple, one-off problems
• Production with fixed solver choice

Modeling Language Trade-offs#

Expressiveness vs Performance#

Algebraic modeling (Pyomo, CVXPY, PuLP):

• More expressive: sum(x[i] for i in range(n))
• Slight overhead: Expression tree construction

Direct API (gurobipy, docplex):

• Less expressive: Must loop manually
• Minimal overhead: Direct solver calls

Overhead magnitude: 5-10% for large models (>10s solve time)

• For small models (<1s), overhead can dominate
• For large models, solver time dominates

Solver Flexibility#

Pyomo example: Same model, different solvers

# Try multiple solvers on same model
for solver_name in ['glpk', 'cbc', 'scip', 'gurobi']:
    solver = SolverFactory(solver_name)
    result = solver.solve(model)
    print(f"{solver_name}: {value(model.obj)}")

Direct API: Requires rewriting for each solver.

Model Maintenance#

Algebraic: Easy to modify structure

# Add new constraint to existing model
model.new_con = Constraint(expr=model.x + model.y >= 5)

Direct API: More programmatic

# Gurobi: add constraint
model.addConstr(x + y >= 5)
model.update()  # Some solvers require explicit update

Hybrid Approaches#

Persistent Solvers (Pyomo)#

Best of both worlds for incremental modeling:

# Create model with Pyomo
model = ConcreteModel()
# ... define model ...

# Create persistent solver instance
solver = SolverFactory('gurobi', solver_io='python')
solver.set_instance(model)

# Solve
solver.solve()

# Modify model
model.new_con = Constraint(expr=model.x <= 10)

# Incremental resolve (fast!)
solver.solve()  # Only communicates changes

CVXPY Parameters#

Separate structure from data for fast re-solves:

# Define parameter (data that changes)
lambda_param = cp.Parameter(nonneg=True)
x = cp.Variable(n)

# Define problem once
objective = cp.Minimize(cp.sum_squares(A @ x - b) + lambda_param * cp.norm(x, 1))
problem = cp.Problem(objective, constraints)

# Solve for different lambda values (fast!)
for lam in [0.01, 0.1, 1.0]:
    lambda_param.value = lam
    problem.solve(warm_start=True)

Comparison Matrix#

Recommendations#

Use Algebraic Modeling When:#

• Model complexity justifies abstraction
• Solver experimentation needed
• Model will be reused/maintained
• Teaching or documentation important
• Not performance-critical

Use Direct API When:#

• Simple, one-off problem
• Maximum performance critical
• Solver-specific features needed
• Incremental model building
• Already committed to specific solver

Use Hybrid (Persistent/Parameters) When:#

• Need both flexibility and performance
• Many similar solves with parameter changes
• Production system with complex model

Comparison with OR-Tools#

OR-Tools uses module-specific APIs (different for LP, CP, routing):

Advantage: Specialized features for each problem type

Disadvantage: Less unified than pure algebraic modeling

Assessment: Good middle ground—more structured than raw API, more specialized than general algebraic modeling.

Migration Paths#

Prototype → Production#

1. Prototype: Algebraic modeling (Pyomo, CVXPY, PuLP)
2. Profile: Identify bottlenecks
3. Optimize:
   • Stay with algebraic if solver time dominates (common)
   • Switch to persistent solvers if model building slow
   • Switch to direct API only if profiling shows clear benefit

Open-Source → Commercial#

Pyomo shines here: Same model code, just change solver

# Development with open-source
solver = SolverFactory('cbc')

# Production with commercial (same model!)
solver = SolverFactory('gurobi')

Key Takeaways#

1. Algebraic modeling overhead small: 5-10% for large problems (solver dominates)
2. Solver flexibility valuable: Pyomo’s 20+ solvers justify abstraction
3. Direct API for performance critical: If profiling shows modeling overhead
4. Hybrid approaches best of both: Persistent solvers, parameters
5. Don’t optimize prematurely: Start with algebraic, profile before switching
6. Problem size matters: Overhead negligible for >10s solves, matters for <1s solves

Recommendation: Default to algebraic modeling (Pyomo, CVXPY, PuLP) unless profiling proves otherwise.

S2-Comprehensive: Solver Performance Comparison#

Overview#

Solver performance varies dramatically based on problem type, size, and structure. This document summarizes benchmark findings and performance characteristics.

Key Finding: Commercial Withdrawal from Public Benchmarks#

Major development (2024): Gurobi (August 2024) and MindOpt (December 2024) withdrew from Mittelmann benchmarks, reducing transparency in commercial vs open-source comparison.

Mittelmann Benchmarks#

Source: Hans Mittelmann, Arizona State University (plato.asu.edu/bench.html)

Coverage:

• Mixed-integer linear programming (MILP)
• Linear programming (LP)
• Network optimization
• Nonlinear programming (NLP)
• Semidefinite programming (SDP)

Interactive visualizations: mattmilten.github.io/mittelmann-plots (shows pairwise time comparisons)

Status (2025): Gurobi and MindOpt results removed. Remaining solvers: SCIP, HiGHS, CBC, GLPK, Xpress, others.

MILP Performance Tiers (Pre-Withdrawal Era)#

Tier 1: Premium Commercial#

Before withdrawal, Gurobi and CPLEX were:

• 5-10x faster than open-source on large instances
• Most consistent across problem classes
• Cost: $10k-100k+ annually (enterprise)
• Academic: Free

When justified:

• Large-scale instances (100k+ variables, 10k+ integer)
• Time-critical applications (real-time optimization)
• Mission-critical production systems

Tier 2: Open-Source Leaders#

SCIP (now Apache 2.0):

• Fastest open-source MILP solver
• Competitive with commercial on many problems
• Active development

HiGHS:

• Rising star (adopted by SciPy 1.6.0+, MATLAB)
• Excellent LP performance
• Good MILP performance
• MIT license, no dependencies

CBC (COIN-OR):

• Mature, stable
• Bundled with PuLP
• 2-5x slower than SCIP on average
• Good for small-medium problems

Tier 3: Basic Open-Source#

GLPK:

• Widely available
• Slower than CBC/SCIP/HiGHS
• Suitable for small problems (<10k variables)

LP Performance#

HiGHS: Competitive with commercial solvers

• Chosen by SciPy and MATLAB as default LP solver
• Signal of high quality

GLOP (Google, in OR-Tools): Good performance for LP

Commercial (Gurobi, CPLEX): Still faster on very large instances, but gap narrowed significantly with HiGHS.

NLP Performance#

Open-Source#

IPOPT (Interior Point OPTimizer):

• Standard open-source NLP solver
• Large-scale capable (10k+ variables)
• Requires linear algebra libraries (HSL, MUMPS, Pardiso)
• Finds local optima (non-convex problems)

Commercial#

Knitro: Premium NLP solver (Artelys)

• Multiple algorithms (interior point, SQP, active set)
• Faster than IPOPT on many problems
• $5k-50k+ annually

SNOPT: Academic/commercial

• Sparse NLP
• Used in aerospace, engineering

Performance gap: Commercial 2-10x faster than IPOPT for large-scale NLP, but IPOPT excellent for open-source.

Convex Optimization Performance#

SOCP/SDP#

MOSEK: Commercial, best performance for conic problems

• Free for academic use
• Industry standard for SOCP/SDP

Open-source:

• SCS: Robust, medium accuracy
• ECOS: Good for small-medium SOCP
• Clarabel: Modern (Rust), becoming default in CVXPY

Performance: MOSEK 2-5x faster, but open-source sufficient for many applications.

Constraint Programming#

OR-Tools CP-SAT: Award-winning

• MiniZinc Challenge medals (2018, 2019, 2020, 2021)
• Competitive with commercial CP solvers
• Free, open-source (Apache 2.0)

Commercial CP:

• IBM CPLEX CP Optimizer
• Gecode

Performance: CP-SAT competitive, sometimes superior.

Problem Size Scaling#

LP#

• Small (<1k variables): All solvers fast (<1 second)
• Medium (1k-100k): HiGHS, commercial excel
• Large (100k-1M+): Interior point methods (HiGHS, Gurobi, CPLEX) scale better than simplex

MILP#

• Small (<100 integer vars): All solvers acceptable
• Medium (100-10k integer): SCIP, HiGHS, commercial
• Large (10k+ integer): Commercial formerly dominated (before withdrawal); SCIP best open-source

NLP#

• Small (<100 vars): scipy.optimize sufficient
• Medium (100-10k): IPOPT, commercial
• Large (10k+): IPOPT, Knitro, SNOPT

Key insight: Problem size alone doesn’t determine hardness. Structure matters more (e.g., network flow LP with 1M variables easier than general MILP with 1k integer variables).

Formulation Impact on Performance#

Example: Same problem, different formulations can solve 10-1000x faster.

Tight formulation: LP relaxation close to integer hull

• Fewer branch-and-bound nodes
• Faster solve

Weak formulation: LP relaxation loose

• Many nodes explored
• Slow solve

Techniques to improve formulation:

• Disaggregation (break large constraints)
• Symmetry breaking (eliminate symmetric solutions)
• Valid inequalities (tighten relaxation)
• Alternative variable representations

Impact often exceeds solver choice: Well-formulated problem on CBC can solve faster than poorly-formulated on Gurobi.

Practical Performance Factors#

1. Presolve#

Modern solvers spend significant time in presolve:

• Can eliminate 50-90% of variables/constraints
• Sometimes solves problem entirely
• Usually worth the time

2. Warm Starting#

MILP: Providing known feasible solution speeds up solving NLP: Good initial guess critical for non-convex

3. Parallelization#

MILP: Scales well to 4-8 cores, diminishing returns beyond 16 NLP: Limited parallelization (inherently sequential)

4. Tolerances#

Relaxing optimality/feasibility tolerances (e.g., 1% vs 0.01%) can yield 10x speedup.

5. Time Limits#

For large MILP, setting time limit and accepting near-optimal solution often practical.

Benchmark Caveats#

1. Problem-dependent: Solver rankings vary by problem class
2. Version-dependent: Solvers improve continuously
3. Tuning: Default parameters may not be optimal
4. Hardware: Benchmark hardware may differ from yours
5. Withdrawn solvers: Gurobi/CPLEX performance no longer publicly tracked

Python Library Overhead#

Modeling language overhead typically <10% of solve time for large problems.

Direct API vs modeling language:

• gurobipy (direct) slightly faster than Pyomo+Gurobi
• Matters for tiny problems (<1s solve) or many sequential solves
• For large problems (>10s), solver dominates

Recommendation: Use algebraic modeling (Pyomo, CVXPY, PuLP) unless profiling shows bottleneck.

Performance Recommendations by Problem Type#

When Commercial Solvers Worth the Cost#

1. Very large MILP: 100k+ variables, 10k+ integer
2. Time-critical: Real-time optimization, tight deadlines
3. ROI justifies: If optimization saves $100k+, $10k solver cost trivial
4. Support needed: Commercial provides professional support

When Open-Source Sufficient#

1. Small-medium problems: Open-source competitive
2. Budget constraints: Free vs $10k-100k
3. Academic/research: Test multiple solvers
4. HiGHS/SCIP for LP/MILP: Closing performance gap

Testing Methodology#

To compare solvers on YOUR problems:

import time
from pyomo.environ import *

def solve_with_solver(model, solver_name):
    solver = SolverFactory(solver_name)
    start = time.time()
    result = solver.solve(model, tee=False)
    elapsed = time.time() - start
    return elapsed, result

# Test multiple solvers
for solver in ['glpk', 'cbc', 'scip', 'gurobi']:
    elapsed, result = solve_with_solver(model, solver)
    print(f"{solver}: {elapsed:.2f}s, obj={value(model.obj)}")