I want to optimize an industrial batch process using pyomo in Python. The objective function is to minimize the cost of the system. In my simplified example, I specify production costs with random values. The uniqueness of the batch process is that once the system is turned on, it must run for four or a multiple of four time steps (multiple batches in a row). The down_time does not have to be a multiple of four. (Only in the symplified problem, actually the downtime should not be longer than two or at least 20 timesteps).
I tried to define a block rule, which should ensures a runtime %4, but it doesn´t work. Still runtimes > 4 appear. If I switch the constraint from >= 4 to == 4 (see below), the optimization doesn´t work at all.
Update:
Everything discussed before is working fine now.
Code:
import pyomo.environ as pyo
import numpy as np
import random
# Define the model
model = pyo.ConcreteModel()
# 1. Define the model, sets, and parameters.
T = 56 # number of timesteps
production_goal = 32
maxPower = 2
minPower = 0
model.timesteps = pyo.RangeSet(1, T)
# Production cost function
production_cost = (np.sin(np.linspace(0, 8*np.pi, T)) + 1) * 50
production_cost = [_ for _ in range(T)]
production_cost = [random.randint(1,100) for _ in range(T)]
# 2. Define the decision variables.
model.production = pyo.Var(model.timesteps, within=pyo.Integers) # production decision
model.running = pyo.Var(model.timesteps, within=pyo.Binary) # plant status
model.start = pyo.Var(model.timesteps, within=pyo.Binary) # start of the plant
# 3. Define the objective function.
def obj_rule(model):
return sum(production_cost[t-1] * model.production[t] for t in model.timesteps)
model.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize)
# 4. Define the constraints.
# Total production constraint
def total_production_rule(model):
return sum(model.production[t] for t in model.timesteps) == production_goal
model.total_production_con = pyo.Constraint(rule=total_production_rule)
# Production is only possible if plant is on
def production_on_rule_1(model, t):
return model.production[t] <= model.running[t] * maxPower
model.production_on_con_1 = pyo.Constraint(model.timesteps, rule=production_on_rule_1)
def production_on_rule_2(model, t):
return model.production[t] >= model.running[t]
model.production_on_con_2 = pyo.Constraint(model.timesteps, rule=production_on_rule_2)
# Running if started in current or preceding 3 periods
def running_rule(model, t):
return model.running[t] <= sum(model.start[tt] for tt in range(t-3, t+1) if tt > 0)
model.running_con = pyo.Constraint(model.timesteps, rule=running_rule)
# Start if no start occurred in preceding 3 periods
def start_rule(model, t):
return model.start[t] <= 1 - sum(model.start[tt] for tt in range(t-3, t) if tt > 0)
model.start_con = pyo.Constraint(model.timesteps, rule=start_rule)
# Shut-down window constraint
M_starts = 5 # max possible starts in 20 periods
def shutdown_window_rule(model, t):
if t > 3: # Limiting timesteps to avoid out-of-bounds errors
return sum(model.start[tt] for tt in range(t+3, t+21) if tt <= T) <= (1 - (model.running[t-1] - model.running[t])) * M_starts
return pyo.Constraint.Skip
model.shutdown_window_con = pyo.Constraint(model.timesteps, rule=shutdown_window_rule)
# Ensure that when the plant starts, it runs for the current and next three periods.
def block_start_rule(model, t):
if t <= T-3:
return model.running[t] + model.running[t+1] + model.running[t+2] + model.running[t+3]>= 4 * model.start[t]
return pyo.Constraint.Skip
model.block_start_con = pyo.Constraint(model.timesteps, rule=block_start_rule)
# 5. Solve the model using Gurobi.
solver = pyo.SolverFactory('gurobi')
results = solver.solve(model)
import matplotlib.pyplot as plt
# 6. Plot the results.
x_values = [model.production[t]() for t in model.timesteps]
y_values = [model.running[t]() for t in model.timesteps]
batch_over_values = [model.start[t]() for t in model.timesteps]
# Create a plot with 3 subplots
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(15, 6))
# Plot the production decision
ax1.step(range(1, T+1), x_values, where='post', label="Production (x)", color='blue')
ax1.set_xticks(range(1, T+1))
ax1.set_ylabel("Value")
ax1.legend()
ax1.grid(True)
# Plot the plant status
ax2.step(range(1, T+1), y_values, where='mid', label="Plant running", color='red', linestyle='--')
ax2.step(range(1, T+1), batch_over_values, where='mid', label="Plant started", color='blue', linestyle=':')
ax2.set_ylabel("Value")
ax2.legend()
ax2.grid(True)
# Plot the production cost
ax3.step(range(1, T+1), production_cost, where='mid', label="Production Cost", color='green', linestyle=':')
ax3.set_xlabel("Timesteps")
ax3.set_ylabel("Value")
ax3.legend()
ax3.grid(True)
# Add a title
plt.suptitle("Steel Plant Production Schedule")
plt.show()
The first thing that I would consider is just making your time periods 4x what they are now and then the whole problem with running for 4 periods goes away, but it appears that you have other constraints that motivate the smaller time step. This leads to some modestly complicated logic to control the running of the plant (as you are discovering) and the constraints to convey that start to grow quickly, but as long as you aren't modeling too large of a timespan, it should be doable.
I'm not sure exactly what x and y are in your model, but you are setting them equal to each other in a constraint, so you really have just 1 variable there. Use clearer names!
Here's a framework idea in pseudocode:
Have a binary variable for when the plant starts and another for if it is running. I have posted several other answers kinda similar to this, so you can search some of my answers for "running" and see what you get. Like this one.
So:
running[t] ∈ {0, 1}
start[t] ∈ {0, 1}
And in order to be running in any period, you need to have started in the current or 3 prior periods.... sounds like a constraint! We can only be running if we started in the preceding 3 periods, or the current one:
running[t] <= sum(start[tt] for tt in range(t-3, t+1) if tt > 0)
for all t
And we cannot arbitrarily "start up" while running, so we need to limit starting after a start occurs. We can only start if no start has occurred in the preceding 3 periods:
start[t] <= 1 - sum(start[tt] for tt in range(t-3, t) if tt > 0)
for all t
And then you can think about the shut-down window. If the system turns off as would be indicated by running[t-1] - running[t] ≥ 0 then we have an exclusion window for new starts from [t+2, t+20]. I thin we can encode that...
sum(start[tt] for tt in range(t+2, t+21)) <= 1 - (running[t-1] - running[t])
for all t ∈ {4, ...}
Try to refactor and incorporate something like that. Keep the model tiny like 9 time periods, print it off and verify your constraints by hand to make sure they are correct. Comment back if you are stuck.
change your constraint to this:
M_starts = 5 # max possible starts in 20 periods
# Shut-down window constraint
def shutdown_window_rule(model, t):
if t > 3: # Limiting timesteps to avoid out-of-bounds errors
return sum(model.start[tt] for tt in range(t+3, t+21) if tt <= T) <= (1 - (model.running[t-1] - model.running[t])) * M_starts
return pyo.Constraint.Skip
What was happening? The earlier version was inadvertently limiting the number of starts in the 20-period window to just one, which was not intended. The constraint needs to be changed to a "big M" constraint so that it basically says "if you are able to start, you can keep restarting in the period...."
Also, the enforcement period was a bit off, it should start at period t+3.
And finally, it needs to cover the latter periods, so the scanning of tt needs to be within the constraint as shown.
This seems to work as desired....
Your next step is (likely) to change your production variable to be a real number and use a big-M constraint to limit that to some max value during the period based on whether the plant is running or not. Good Luck!