Optimizing Production Schedules with Mixed-Integer Linear Programming

Imagine you run a manufacturing facility that produces widgets 24/7. Your electricity costs vary by the hour — cheap at night, expensive during peak demand. Your machine can run at different production rates, but higher speeds consume more power (and not always proportionally). You have a daily production target to meet.

The question: When should you run the machine, and at what rate, to minimise your energy costs while meeting your production goals?

This is a classic day-ahead scheduling problem, and it’s exactly the kind of challenge that Mixed-Integer Linear Programming (MILP) excels at solving. In this blog, I’ll walk through building a complete MILP optimiser in Python using Pyomo and the CBC solver. Along the way, we’ll discover some fascinating insights about how equipment efficiency curves fundamentally shape optimal production strategies. You can see the code here in GitHub.

What is MILP?

Mixed-Integer Linear Programming combines two types of decision variables:

• Integer variables: Discrete choices (e.g., machine ON/OFF, which shift to schedule)
• Continuous variables: Quantities that can vary smoothly (e.g., production rate, power consumption)

Unlike pure Linear Programming (LP), MILP can model real-world constraints like binary on/off decisions, minimum up/down times (can’t cycle too quickly), startup costs, and discrete operating modes. The “Linear” part means the objective and constraints must be linear — but we’ll see how to handle nonlinear power curves using clever approximations.

The problem: production scheduling with nonlinear power consumption

Inputs

• 24-hour electricity prices (€/kWh) — varying from €0.09 to €0.22
• Production target: 200 widgets (must produce exactly this amount)
• Machine constraints: rate limits 10–100 widgets/hr when ON; minimum up time 2 hours; minimum down time 1 hour

The power consumption curve is nonlinear:

power(kW) = a × rate² + b × rate + c

where a is the quadratic coefficient (inefficiency at high rates), b the linear coefficient (direct energy per widget), and c the base load (fixed overhead when ON).

Output

An optimal 24-hour schedule specifying which hours to run the machine (ON/OFF), the production rate for each hour, and the total energy cost.

The MILP formulation

Decision variables

Binary (integer) variables — 72 total (24 hours × 3):

m.on[t]    # 1 if machine ON at hour t, 0 otherwise
m.start[t] # 1 if machine starts at hour t
m.stop[t]  # 1 if machine stops at hour t

Continuous variables — 48 total (24 hours × 2):

m.rate[t]  # Production rate (widgets/hour) at hour t
m.power[t] # Power consumption (kW) at hour t

Key constraints

1. Link rate to on/off state:

rate[t] >= rate_min * on[t]  # When ON: rate >= 10
rate[t] <= rate_max * on[t]  # When OFF: rate = 0

2. Startup/shutdown logic:

start[t] >= on[t] - on[t-1]  # Detects OFF→ON transitions
stop[t]  >= on[t-1] - on[t]  # Detects ON→OFF transitions

3. Minimum up/down time:

# If we start at t, must stay on for min_up_hours
sum(on[k] for k in range(t, t+min_up_hours)) >= min_up_hours * start[t]

# If we stop at t, must stay off for min_down_hours
sum(1-on[k] for k in range(t, t+min_down_hours)) >= min_down_hours * stop[t]

4. Daily production target:

sum(rate[t] for t in 0..23) == 200  # Equality mode

5. Piecewise-linear power curve (the clever part): since power = a×rate² + b×rate + c is nonlinear, we approximate it using piecewise-linear segments with SOS2 (Special Ordered Set type 2) variables:

m.pw = pyo.Piecewise(
    m.T,                    # Index: hours 0-23
    m.power,                # Dependent variable: power
    m.rate,                 # Independent variable: rate
    pw_pts=breakpoints,     # e.g., [0, 9, 18, 27, ..., 100]
    f_rule=power_values,    # Power at each breakpoint
    pw_constr_type="EQ",    # Equality: power = f(rate)
    pw_repn="SOS2",         # Use SOS2 for MILP efficiency
)

This creates 12 piecewise-linear segments that closely approximate the quadratic curve while keeping everything linear for the MILP solver.

Objective function

minimize: energy_cost + startup_penalties

where:
  energy_cost       = sum(price[t] * power[t] for t in 0..23)
  startup_penalties = startup_cost * sum(start[t] for t in 0..23)

The fascinating role of the quadratic coefficient

By varying just the quadratic coefficient a in the power curve, we can simulate different types of machinery — and the optimiser adapts its strategy dramatically.

Experiment 1: high quadratic penalty (a = 0.015)

Power curve: power = 0.015×rate² + 0.6×rate + 5.0. At maximum rate (100 widgets/hr) the quadratic term is 150 kW (70% of total), linear 60 kW, base 5 kW — total 215 kW. This represents equipment with severe inefficiencies at high speeds (e.g. pumps with quadratic drag losses).

Hours ON: 0-7, 15-20, 23 (15 hours total)
Max rate: 27 widgets/hr (conservative)
Total cost: €35.62

The optimiser spreads production across many hours at low rates to avoid the crushing quadratic penalty.

Experiment 2: low quadratic penalty (a = 0.001)

Power curve: power = 0.001×rate² + 0.6×rate + 5.0. At maximum rate the total is just 75 kW. This represents highly scalable equipment with mostly linear power consumption.

Hours ON: 2-6 (5 hours only!)
Max rate: 60 widgets/hr (aggressive burst)
Total cost: €29.25 (18% savings!)

The optimiser concentrates all production in the cheapest hours (4–5am at €0.09/kWh) and runs at maximum feasible rates: a “sprint during cheap hours, stop during expensive hours” strategy.

Experiment 3: medium quadratic penalty (a = 0.005)

Power curve: power = 0.005×rate² + 0.6×rate + 5.0 — total 115 kW at max rate.

Hours ON: 0-6, 18-19 (9 hours)
Max rate: 36 widgets/hr (balanced)
Total cost: €31.75

A perfect middle ground — moderate rates, focused on cheaper hours.

The strategy spectrum

Real-world solver performance

Using the CBC (COIN-OR Branch and Cut) solver:

Problem size: 191 constraints, 383 variables (71 binary)
LP relaxation: €29.25 (lower bound)
Cutting planes: 13 cuts added (Gomory, Probing, MIR)
Integer solution: €29.25 (found in 0 nodes!)
Solve time: 0.01 seconds

Key insight: the cutting planes were so effective that the LP relaxation bound matched the integer solution exactly — no branch-and-bound exploration needed. This is why MILP solvers are so powerful for well-structured problems.

Practical takeaways

1. Equipment efficiency curves matter — a lot

The quadratic coefficient a acts as a “burst penalty” dial: low a (scalable equipment) favours a burst strategy during cheap periods; high a favours spreading across more hours. Before optimising, measure your actual equipment’s power curve:

curve:
  type: breakpoints
  rate_points:  [0, 10, 20, 40, 60,  80,  100]
  power_points: [0, 12, 25, 52, 85, 130, 195]  # From real measurements

2. MILP can handle complexity

This problem has binary on/off decisions, minimum up/down constraints, nonlinear power curves (approximated), and hourly-varying constraints and prices — yet it solves in 0.01 seconds on a laptop.

3. The base load (c) creates trade-offs

The fixed 5 kW base load when ON makes short, low-rate runs inefficient (overhead dominates), encouraging either running longer or at higher rates — the “sprint vs. marathon” dynamic.

4. Constraints shape solutions in non-obvious ways

The min_up_hours=2 constraint forces at least 2-hour runs, preventing rapid on/off cycling that would otherwise chase every price fluctuation (unrealistic for real equipment).

The code

The complete optimiser is ~400 lines of Python and supports polynomial or breakpoint power curves, CSV input for prices and hourly constraints, configurable min up/down times, ramp rates and startup costs, and both equality and minimum production targets. Key dependencies and run command:

pip install pyomo pandas pyyaml
sudo apt install coinor-cbc   # or: brew install cbc

python producer_milp.py \
  --config config.yaml \
  --prices prices.csv \
  --out schedule.csv

Conclusion

By modelling our manufacturing problem as a MILP, we found 18% cost savings by tuning equipment strategy to match its efficiency characteristics, discovered how quadratic efficiency curves change optimal scheduling, and solved a complex 24-hour problem in 10 milliseconds. The real magic happens when you combine domain knowledge (equipment physics, operational constraints), mathematical modelling (the MILP formulation), and modern solvers (CBC, Gurobi, CPLEX). Whether you’re scheduling production, optimising energy systems, planning logistics, or routing vehicles — MILP is likely the right tool for the job.

The full code is available on GitHub under the MIT License. Copyright © 2026 FullStackEnergy.com.