gurobi_optimizer.py

"""
Solve example problem as MIP using Gurobi.
"""
import numpy as np
import pandas as pd
import gurobipy
from gurobipy import GRB
from interface import ProductionProblem


class GurobiSolver:
    """MIP solver using Gurobi.

    Parameters
    ----------
    time_limit : float
        Max allowed runtime for the solver in seconds
    """

    def __init__(self, time_limit=60):
        self.model = gurobipy.Model("hcem-mips")
        self.model.Params.TimeLimit = time_limit

    def calc_plan(self, problem: ProductionProblem):
        """Interface compatible with the one used in evaluation.py

        I'm not a fan, it would be better to pass `problem` on init,
        and to create object state variables there instead of like
        here in an instance method.
        """
        self.problem = problem
        self.make_gurobi_model()
        self.model.optimize()

        if True:
            # Convert to format that the Evaluator needs
            return self.solution.rename(
                columns={"milling": "on", "time": "datetime"}
            ).set_index("datetime").assign(mill_id=1, cement_type_id=1, silo_id=1)

    def make_gurobi_model(self):
        """Translate problem into Gurobi model object"""
        time_interval = 1
        price = np.array(self.problem.price_forecast)
        n_time = len(price)
        final_stock = self.problem.target_fill * self.problem.capacity
        consumption = np.array(self.problem.demand_forecast)

        # Variable: is machine on in time interval?
        mill_running = self.model.addMVar(n_time, vtype=GRB.BINARY, name="mill_running")

        # Variable: silo fill level
        silo_level = self.model.addMVar(
            n_time,
            lb=self.problem.capacity * self.problem.min_fill,
            ub=self.problem.capacity,
            vtype=GRB.CONTINUOUS,
            name="f",
        )

        # Objective: minimize electricity price
        self.model.setObjective(
            mill_running @ price * time_interval * self.problem.power,
            sense=GRB.MINIMIZE,
        )

        # Contraint on initial and final stock
        self.model.addConstr(silo_level[0] == self.problem.initial_stock)
        self.model.addConstr(silo_level[-1] >= final_stock)

        # Constraint to declare how stock, production, consumption are related
        for i in range(n_time - 1):
            production = mill_running[i] * time_interval * self.problem.throughput
            self.model.addConstr(
                silo_level[i + 1] == silo_level[i] + production - consumption[i]
            )

        # Minimum mill runtime constraint
        # TODO: meditate and debug this part.
        # At the moment just copy & paste & adapt from mip_optimizer.py
        min_runtime = int(self.problem.min_runtime)
        for i in range(0, n_time):
            for k in range(1, min_runtime):
                if i + k < n_time:
                    if i == 0:
                        self.model.addConstr(mill_running[i] <= mill_running[i+k])
                    else:
                        self.model.addConstr((1-mill_running[i-1])+mill_running[i]-1 <= mill_running[i+k])
                else:
                    # not sure about this condition (edge case at the end of schedule)
                    self.model.addConstr(mill_running[i] <= mill_running[i-k])

        # TBD: there doesn't seem to be a way to retrieve an MVar yet from the model in Gurobi
        # https://support.gurobi.com/hc/en-us/community/posts/360061782092-Retrieving-MVar-from-Model-read-from-file
        # https://www.gurobi.com/documentation/9.0/quickstart_mac/py_reporting_results_attri2.html
        # As a workaround, we store the MVar objects here
        # so that we can access them later when evaluating the solution
        self.model_variables = {
            "mill_running": mill_running,
            "silo_level": silo_level,
        }

    @property
    def solution(self):
        """Translate Gurobi model into Python solution data structure"""
        milling = self.model_variables["mill_running"].X.astype(bool)
        stock = self.model_variables["silo_level"].X
        consumption = np.array(self.problem.demand_forecast)
        production = milling * self.problem.throughput
        # Cross-check the stock constraint implementation above - should be identical
        stock2 = self.problem.initial_stock + production.cumsum() - consumption.cumsum()

        return pd.DataFrame(
            {
                "time": self.problem.timestamp,
                "price": self.problem.price_forecast,
                "consumption": consumption,
                "milling": milling,
                "stock": stock,
                "stock2": stock2,
            }
        )

    def plot_solution(self):
        import matplotlib.pyplot as plt

        df = self.solution

        fig, axes = plt.subplots(3, 1, sharex=True, figsize=(12, 10))

        ax = axes[0]
        ax.plot(df.time, df.price, label="not milling")
        ax.plot(df.time[df.milling], df.price[df.milling], "ro", ms=3, label="milling")
        ax.set_ylabel("price (Euro)")
        ax.grid()
        ax.legend()

        ax = axes[1]
        ax.plot(df.time, df.stock, label="stock")
        ax.plot(df.time, df.stock2, ls="dotted", label="stock2")
        ax.plot(df.time[df.milling], df.stock[df.milling], "ro", ms=3)
        ax.axhline(self.problem.min_fill * self.problem.capacity, label="min allowed")
        ax.axhline(self.problem.capacity, label="max allowed")
        ax.set_ylim(0, None)
        ax.set_xlabel("Date")
        ax.set_ylabel("stock (tonnes)")
        ax.grid()
        ax.legend()

        ax = axes[2]
        ax.plot(df.time, df.consumption)
        ax.set_ylabel("consumption (tonnes)")

        fig.savefig("gurobi_optimizer.png")


def main():
    problem = ProductionProblem.read(
        params="params/params.json",
        forecast="input_data/geseke_elsa_eh1_cd_day_start.csv",
    )

    solver = GurobiSolver()
    solution = solver.calc_plan(problem)
    print(solution)
    # print(solver.model)
    # print(solver.solution)

    solver.plot_solution()


if __name__ == "__main__":
    main()