Solver Max - Schedule staff with enumerated shifts, Pyomo

3 April 2023

Schedule enumerated

In the previous article, we built a staff scheduling optimization model in Excel, using OpenSolver. That model enumerates all possible shift patterns and then decides how many of each shift pattern to use so that we meet the various constraints at least cost.

Our focus in this article is on replicating the Excel model by translating it into Python, using the Pyomo library. That involves formalizing the mathematical formulation, creating appropriate data structures, and representing the formulation in Pyomo syntax. Both models use the CBC solver, so the differences between the models aren't in the solving process, but rather in the build process.

Excel can be a great tool for some types of optimization modelling, either for prototyping or as a production tool. Python offers a wider range of optimization options, and it is usually better for models that have varying data sizes or need to operate at a large scale. The advantages and disadvantages of modelling in Excel and Python are discussed in more detail in our article Optimization in Excel vs Python.

Download the model

The model is available on GitHub.

Situation

The modelling situation is identical to that in the previous article. That is, our business opens at 6:00am and closes at 6:00pm. We need to create a shift schedule for our staff to meet the demand, expressed as the minimum number of staff required in each half hour. Demand is low at the start of the working day, rises to a peak at around noon, then declines towards the end of the working day. Our objective is to minimize the total wage cost for the day.

We have three categories of staff. Each category may be assigned to shifts that start at the beginning of any half-hour, provided that:

Full-time staff work 8 hours per day, plus a 1 hour unpaid break in the middle of their shift. Up to 7 of these staff are available, at a wage of $22.00/hour each.
Some part-time staff work 6 hours per day, plus a 0.5 hour unpaid break in the middle of their shift. Up to 6 of these staff are available, at a wage of $21.00/hour each.
Some part-time staff work 4 hours per day, with no unpaid break. Up to 14 of these staff are available, at a wage of $20.00/hour each.

Model design

Get the data

Our Python model needs some data. We already have all the data we need in the Excel version of this model, so we'll use that file for data input. A convenient way to read data from Excel is to use named ranges – our Excel file already has all the data ranges named, so we can use them to access the data.

Define the sets

Excel and Python have very different ways of representing an optimization model. Excel uses a grid structure, which implicitly defines data dimensions. In Pyomo, we need to explicitly define a Set for each dimension.

Examining the Excel workbook, we see that there are three implicit sets used to define the data. That is:

People. The categories of people who work in our business, defined in terms of the hours they work: 8h, 6h, or 4h per day.
Shifts. An enumerated list of shifts that our people can work. In this model, there are 36 possible shifts, representing various start times for each of the people categories.
Slots. We need to meet demand for specific time slots over the working day, beginning at 6:00am through to 5:30pm in half-hour steps.

We will need to define each of these sets in our Pyomo model. Then we'll need to translate the raw data into Pyomo objects, mostly Param parameters, defined over these sets.

Define the variables

In the Solver/OpenSolver dialog, we can see that there is one set of variables, vAllocation, defined over the Shifts dimension. These variables represent the number of people allocated to each of the enumerated shifts.

Define the constraints

The extended Excel model has five constraints, as shown in the Solver or OpenSolver dialog. Similarly, we'll need to define five constraints in Pyomo, representing:

Demand. Ensure we have the required number of staff, in total, in each time slot.
Availability. Number of staff of each category is within their availability.
Minimum people. Number of staff of each category is at least the minimum required.
Shift must use. Must use specific shifts, if required.
Surplus upper bound. The number of surplus staff in each time slot is no more than the maximum allowed.

In the Excel model, we construct each constraint using a combination for formulae and selections in the Solver/OpenSolver dialog. In Pyomo, we'll need to represent each constraint using Python code.

Define the objective function

In the Solver/OpenSolver dialog, the objective function refers to a single cell. The value in that cell is the result of a series of calculations that involve many SUM and SUMIFS functions, along with other calculations. We'll need to untangle those calculations and translate them into Pyomo syntax – a process that, given the different perspectives of Excel and Python models, is often not straightforward.

Model formulation

A key in translating the model is to formalize the mathematical formulation. To some extent we can ignore the maths when building a model in Excel. But having a formal formulation is very useful when writing the Pyomo version, as the Python code is an almost literal translation of the formulation.

The formulation for our model is shown in Figure 1.

Figure 1. Model formulation

\begin{alignat}{1} &\text{Objective} \\ & \quad \text{Minimize} \ TotalCost &= \quad & \sum_{p=1}^m \sum_{s=1}^n \text{ShiftSlots}_{s, t} \times Allocation_{s} \\ &&& \times \text{Wage}_{p} \times \text{TimePeriod}\tag{1} \\ &&& \text{if } \text{ShiftLabel}_{s} = \text{PeopleLabel}_{p} \\ \\ &\text{Subject to} \\ & \quad \sum_{s=1}^n \text{ShiftSlots}_{s, t} \times Allocation_{s} &\ge \quad & \text{SlotRequired}_{t} \quad &\forall \ t \in T \tag{2}\\ & \quad \sum_{s=1}^n Allocation_{s} &\le \quad & \text{Available}_{p} \quad \text{if } \text{ShiftLabel}_{s} = \text{PeopleLabel}_{p} \quad &\forall \ p \in P \tag{3}\\ & \quad \sum_{s=1}^n Allocation_{s} &\ge \quad & \text{PeopleMin}_{p} \quad \text{if } \text{ShiftLabel}_{s} = \text{PeopleLabel}_{p} \quad &\forall \ p \in P \tag{4}\\ & \quad Allocation_{s} &\ge \quad & \text{ShiftRequired}_{s} \quad &\forall \ s \in S \tag{5}\\ & \quad \sum_{s=1}^n \text{ShiftSlots}_{s, t} \times Allocation_{s} \quad &\le \quad & \text{MaxSurplus} + \text{SlotRequired}_{t} \quad &\forall \ t \in T \tag{6}\\ \\ &\text{Variables} \\ & \quad Allocation_{s} &= \quad & \text{Number of people per shift} \quad &\forall \ s \in S \tag{7}\\ \\ &\text{Data} && \\ & \quad \text{TimePeriod} &= \quad &\text{Hours per time slot} \quad \tag{8} \\ & \quad \text{Wage}_{p} &= \quad &\text{Wage rate per person, \$/hour} \quad &\forall \ p \in P \tag{9}\\ & \quad \text{ShiftSlots}_{s, t} \quad &= \quad &\text{Definition of each shift by time slot} \quad &\forall \ s \in S, \forall \ t \in T \tag{10}\\ & \quad \text{ShiftLabel}_{s} &= \quad &\text{Person category of each shift, non-unique} \quad &\forall \ s \in S \tag{11}\\ & \quad \text{PersonLabel}_{p} &= \quad &\text{Label for each person category} \quad &\forall \ p \in P \tag{12}\\ & \quad \text{SlotRequired}_{t} &= \quad &\text{Number of people required in each time slot} \quad &\forall \ t \in T \tag{13}\\ & \quad \text{Available}_{p} &= \quad &\text{Number of each person category available} \quad &\forall \ p \in P \tag{14}\\ & \quad \text{PeopleMin}_{p} &= \quad &\text{Minimum people per category} \quad &\forall \ p \in P \tag{15}\\ & \quad \text{ShiftRequired}_{s} &= \quad &\text{Minimum people allocated to each shift} \quad &\forall \ s \in S \tag{16}\\ & \quad \text{MaxSurplus} &= \quad &\text{Maximum surplus people (all slots)} \quad \tag{17} \\ \\ &\text{Sets} \\ & \quad P &\ &\text{People category} \tag{18} \\ & \quad S &\ &\text{Shifts} \tag{19} \\ & \quad T &\ &\text{Slots} \tag{20} \\ \\ &\text{Dimensions} \\ & \quad m &\ &\text{Number of people categories} \tag{21} \\ & \quad n &\ &\text{Number of shifts} \tag{22} \\ & \quad k &\ &\text{Number of slots} \tag{23} \\ \end{alignat}

Model implementation

Import dependencies

The first task is to import the libraries that are needed for our program, as shown in Figure 2. The dependencies include the openpyxl library, which we'll use to read data from our Excel workbook.

Figure 2. Import dependencies

import pyomo.environ as pyo
import pandas as pd
import os.path
from openpyxl import load_workbook
from openpyxl.utils.cell import range_boundaries

Globals

We have a few global objects that will make things easier, as shown in Figure 3. This includes Model, which is the primary object used for building the model.

As an aside, we acknowledge that programming purists will frown upon the use of global values. However, these values are the inputs that the user is most likely to change, so it is convenient to have them collated near the top of the program and available throughout the program.

Figure 3. Global objects

Model = pyo.ConcreteModel(name = 'Schedule with enumerated shifts')
ExcelFile = 'Schedule-staff-with-enumerated-shifts.xlsx'
Worksheet = 'Model 1'
Engine = 'cbc'   # cbc, glpk, or appsi_highs
TimeLimit = 60   # seconds

Excel workbook generic loader

Since we'll get our data from named ranges in an Excel workbook, it is helpful to have a generic function to perform this task, as shown in Figure 4. We just supply this function an Excel workbook file name, a worksheet name, and a range name – it returns a Pandas dataframe containing that data.

Figure 4. Generic loader from Excel file, given worksheet and named range

def LoadFromExcel(ExcelFile, Worksheet, Range):
    Workbook = load_workbook(filename=ExcelFile, read_only = True)
    Sheet = Workbook[Worksheet]
    NamesRanges = Workbook.defined_names[Range].destinations
    for Title, Coord in NamesRanges:
        MinCol, MinRow, MaxCol, MaxRow = range_boundaries(Coord)
        Data = Sheet.iter_rows(MinRow, MaxRow, MinCol, MaxCol, values_only = True)
    ExtractedData = pd.DataFrame(Data)
    return ExtractedData

Get data

Figure 5 shows the process for getting each data item, given the range name where it is located in the workbook. Each of these dataframes is just a temporary placeholder – we'll put them into the Model object in the next section.

Figure 5. Get data from workbook

def GetData(ExcelFile, Worksheet):
    Available = LoadFromExcel(ExcelFile, Worksheet, 'dAvailable')
    People = LoadFromExcel(ExcelFile, Worksheet, 'dPeople')
    PeopleHourly = LoadFromExcel(ExcelFile, Worksheet, 'dPeopleHourly')
    PeopleMin = LoadFromExcel(ExcelFile, Worksheet, 'dPeopleMin')
    SlotRequired = LoadFromExcel(ExcelFile, Worksheet, 'dPeopleRequired')
    Shifts = LoadFromExcel(ExcelFile, Worksheet, 'dShifts')
    ShiftSlots = LoadFromExcel(ExcelFile, Worksheet, 'dShiftSlots')
    Slots = LoadFromExcel(ExcelFile, Worksheet, 'dSlots')
    SurplusMax = LoadFromExcel(ExcelFile, Worksheet, 'dSurplusMax')
    ShiftRequired = LoadFromExcel(ExcelFile, Worksheet, 'dShiftRequired')
    TimePeriod = LoadFromExcel(ExcelFile, Worksheet, 'dTimePeriod')
    return Available, People, PeopleHourly, PeopleMin, SlotRequired, Shifts, ShiftSlots, Slots, SurplusMax, ShiftRequired, TimePeriod

Declarations

The data declarations are shown in Figure 6. This is where we translate the Excel data into the Model object.

We divide the declarations into three sections:

Single parameter values. This includes things like the maximum surplus assumption (a single number) and the solver name (a string).
Sets. The sets correspond to the indices in the formulation. That is, People, Shifts, and Slots. Note that we're using Pyomo's Set, which is 1-based, rather than Python's 0-based set.
Parameters. In populating the set-based parameters, our general approach is to declare a parameter as mutable, and then separately populate it with data. Note that we need to adjust some of the indices to account for the difference between 0-based dataframes and 1-based Pyomo Sets.

Figure 6. Declarations

def DefineData():
    Available, People, PeopleHourly, PeopleMin, SlotRequired, Shifts, ShiftSlots, Slots, SurplusMax, ShiftRequired, TimePeriod = GetData()

    # Single parameter values
    Model.TimePeriod = pyo.Param(within = pyo.NonNegativeReals, initialize = TimePeriod[0][0].item())   # Hours per time slot
    Model.MaxSurplus = pyo.Param(within = pyo.Any, initialize = SurplusMax[0][0].item())   # Maximum surplus people (all slots)
    Model.Engine = pyo.Param(within = pyo.Any, initialize = Engine)   # Name of solver engine (global)
    Model.TimeLimit = pyo.Param(within = pyo.NonNegativeReals, initialize = TimeLimit)   # Solver time limit (global)

    # Sets
    Model.People = pyo.Set(initialize = range(1, len(People) + 1))   # Categories of people
    Model.Shifts = pyo.Set(initialize = range(1, len(Shifts) + 1))   # Enumerated shifts
    Model.Slots = pyo.Set(initialize = range(1, len(Slots.columns) + 1))   # Time slots

    # Parameters using sets
    Model.PeopleLabel = pyo.Param(Model.People, within = pyo.Any, mutable = True)   # Label for each person category (unique)
    Model.Available = pyo.Param(Model.People, within = pyo.NonNegativeIntegers, mutable = True)   # Number of each person category available
    Model.Wage = pyo.Param(Model.People, within = pyo.NonNegativeReals, mutable = True)   # Wage rate per person, $/hour
    Model.PeopleMin = pyo.Param(Model.People, within = pyo.NonNegativeIntegers, mutable = True)   # Minimum number of each person category to be used
    for p in Model.People:    
        Model.PeopleLabel[p] = People[0][p - 1]   # Single column
        Model.Available[p] = Available[0][p - 1]
        Model.Wage[p] = PeopleHourly[0][p - 1]
        Model.PeopleMin[p] = PeopleMin[0][p - 1]

    Model.ShiftLabel = pyo.Param(Model.Shifts, within = pyo.Any, mutable = True)   # Person category of each shift (non-unique)
    Model.ShiftRequired = pyo.Param(Model.Shifts, within = pyo.NonNegativeIntegers, mutable = True)   # Minimum number of people allocated to each shift
    for s in Model.Shifts:
        Model.ShiftLabel[s] = Shifts[0][s - 1]   # Single column
        Model.ShiftRequired[s] = ShiftRequired[0][s - 1]

    Model.SlotLabel = pyo.Param(Model.Slots, within = pyo.Any, mutable = True)   # Label for each time slot (output only)
    Model.SlotRequired = pyo.Param(Model.Slots, within = pyo.NonNegativeIntegers, mutable = True)   # Number of people required in each time slot
    for t in Model.Slots:
        Model.SlotLabel[t] = Slots[t - 1][0]   # Single row
        Model.SlotRequired[t] = SlotRequired[t - 1][0]

    Model.ShiftSlots = pyo.Param(Model.Shifts, Model.Slots, within = pyo.NonNegativeIntegers, mutable = True)   # Definition of each shift by time slot (0, 1)
    for s in Model.Shifts:
        for t in Model.Slots: 
        Model.ShiftSlots[s, t] = ShiftSlots[t - 1][s - 1] # Translate from 0-base to 1-base

Check the data

It can be useful to check that we've correctly translated the Excel data into the Model object. Checking that the data is as expected can save a lot of frustration when building the rest of the model, so we do that in the function shown in Figure 7. This function has a Boolean parameter that determines whether it outputs anything – effectively allowing us to skip the check.

Figure 7. Check data

def CheckData(Display):
    if Display:
        print('Data check')
        print('----------\n')

        print('Time period: ', pyo.value(Model.TimePeriod))
        print('Max surplus: ', pyo.value(Model.MaxSurplus))
        print('Engine: ', pyo.value(Model.Engine))
        print('TimeLimit: ', pyo.value(Model.TimeLimit))

        print('\nPeople set:', end = ' ')
        for p in Model.People:
            print(p, end = ' ')

        print('\nShifts set:', end = ' ')
        for s in Model.Shifts: 
            print(s, end = ' ')

        print('\nSlots set:', end = ' ')
        for t in Model.Slots: 
            print(t, end = ' ') 

        print('\nPeople labels:', end = ' ')
        for p in Model.People:
            print(pyo.value(Model.PeopleLabel[p]), end = ' ')

        print('\nPeople attributes:')
        for p in Model.People:
            print(f'{pyo.value(Model.PeopleLabel[p])}:  Available {pyo.value(Model.Available[p]):4}  Min {pyo.value(Model.PeopleMin[p]):4}  Wage {pyo.value(Model.Wage[p]):6.2f}')

        print('\nShift labels:', end = ' ')
        for s in Model.Shifts:
            print(pyo.value(Model.ShiftLabel[s]), end = ' ')

        print('\nShift required:', end = ' ')
        for s in Model.Shifts:
            print(pyo.value(Model.ShiftRequired[s]), end = ' ')

        print('\nSlot labels:', end = ' ')
        for t in Model.Slots:
            print(pyo.value(Model.SlotLabel[t]), end = ' ')

        print('\nSlot required:', end = ' ')
        for t in Model.Slots:
            print(pyo.value(Model.SlotRequired[t]), end = ' ')

        print('\nShift slots:')
        for s in Model.Shifts:
            for t in Model.Slots: 
                print(pyo.value(Model.ShiftSlots[s, t]), end = ' ')
            print()

        print()

Define the model

The model definition is shown in Figure 8.

There are five constraint definitions, corresponding to the five constraints in the formulation:

Model.PeopleRequired. Ensure we have the required number of staff, in total, in each time slot.
Model.PeopleAvailable. Number of staff of each category is within their availability.
Model.PeopleLB. Number of staff of each category is at least the minimum required.
Model.ShiftMustUse. Must use specific shifts, if required.
Model.SlotSurplus. The number of surplus staff in each time slot is no more than the maximum allowed. In Equation (6) of the formulation we moved the Model.SlotRequired term to the right-hand-side for presentation reasons.

We provide two alternative objective function definitions. Model.TotalCostAlt is a direct translation of what the spreadsheet formulae do. Model.TotalCost is entirely equivalent, but we think it is easier to understand.

Note that the objective function and two of the constraints include the expression if pyo.value(Model.ShiftLabel[s]) == pyo.value(Model.PeopleLabel[p]) within the sum. This test is used in the formulation to represent summing by the people categories 8h, 6h, and 4h – which we did in the spreadsheet using separate rows.

Figure 8. Define the model

def DefineModel():
    Model.Allocation = pyo.Var(Model.Shifts, domain = pyo.NonNegativeIntegers)   # Number of each shift to use

    def rule_demand(Model, t):   # Ensure we have the required number of staff, in total, in each time slot
        return sum(Model.ShiftSlots[s, t] * Model.Allocation[s] for s in Model.Shifts) >= Model.SlotRequired[t]
    Model.PeopleRequired = pyo.Constraint(Model.Slots, rule = rule_demand)

    def rule_available(Model, p):   # Number of staff of each category is within their availability
        return sum(Model.Allocation[s] for s in Model.Shifts if pyo.value(Model.ShiftLabel[s]) == pyo.value(Model.PeopleLabel[p])) <= Model.Available[p]
    Model.PeopleAvailable = pyo.Constraint(Model.People, rule = rule_available)

    def rule_peoplemin(Model, p):   # Number of staff of each category is at least the minimum required
        return sum(Model.Allocation[s] for s in Model.Shifts if pyo.value(Model.ShiftLabel[s]) == pyo.value(Model.PeopleLabel[p])) >= Model.PeopleMin[p]
    Model.PeopleLB = pyo.Constraint(Model.People, rule = rule_peoplemin)

    def rule_mustuseshift(Model, s):   # Must use specific shifts, if required
        return Model.Allocation[s] >= Model.ShiftRequired[s]
    Model.ShiftMustUse = pyo.Constraint(Model.Shifts, rule = rule_mustuseshift)

    def rule_surplus(Model, t):   # The number of surplus staff in each time slot is no more than the maximum allowed
        return sum(Model.ShiftSlots[s, t] * Model.Allocation[s] for s in Model.Shifts) - Model.SlotRequired[t] <= Model.MaxSurplus
    Model.SlotSurplus = pyo.Constraint(Model.Slots, rule = rule_surplus)

    # def rule_ObjAlt(Model):   # Alternative, equivalent way of writing the objective function
    #     Cost = 0
    #     for p in Model.People:
    #         Cost += sum(sum(Model.ShiftSlots[s, t] * Model.Allocation[s] for s in Model.Shifts if pyo.value(Model.ShiftLabel[s]) == pyo.value(Model.PeopleLabel[p])) for t in Model.Slots) * pyo.value(Model.Wage[p]) * Model.TimePeriod
    #     return Cost
    # Model.TotalCostAlt = pyo.Objective(rule = rule_ObjAlt, sense = pyo.minimize)

    def rule_Obj(Model):   # Minimize total cost of staff
        Cost = 0
        for p in Model.People:
            for s in Model.Shifts:
                if pyo.value(Model.ShiftLabel[s]) == pyo.value(Model.PeopleLabel[p]):
                    Cost += sum(Model.ShiftSlots[s, t] * Model.Allocation[s] for t in Model.Slots) * pyo.value(Model.Wage[p])
        Cost = Cost * Model.TimePeriod
        return Cost
    Model.TotalCost = pyo.Objective(rule = rule_Obj, sense = pyo.minimize)

Solve model

We're now ready to solve the model, as shown in Figure 9. We've allowed for a choice between different solvers, each of which specifies the time limit differently. In addition, when calling this function, we can specify whether we want the solver to produce verbose output as it works on the solution.

Note that, as of March 2023, the Pyomo interface to the HiGHS solver has a bug where mutable parameters with zero values produce an error. Therefore, the HiGHS solver doesn't currently work for this model – though hopefully it will soon.

Figure 9. Solve model

def SolveModel(Verbose):
    Solver = pyo.SolverFactory(pyo.value(Model.Engine))

    if pyo.value(Model.Engine) == 'cbc':
        Solver.options['seconds'] = pyo.value(Model.TimeLimit)
    elif pyo.value(Model.Engine) == 'glpk':
        Solver.options['tmlim'] = pyo.value(Model.TimeLimit)
    elif pyo.value(Model.Engine) == 'appsi_highs':
        Solver.options['time_limit'] = pyo.value(Model.TimeLimit)

    Model.dual = pyo.Suffix(direction = pyo.Suffix.IMPORT)

    Results = Solver.solve(Model, load_solutions = False, tee = Verbose)
    return Results

Process results

After the solve, we process the results, as shown in Figure 10. Note that in the solve process above, we defer loading the results in case there are no results to load (because the solver didn't return a result, perhaps due to an error).

Figure 10. Process results

def ProcessResults():
    WriteSolution = False
    Optimal = False
    LimitStop = False
    Condition = Results.solver.termination_condition

    if Condition == pyo.TerminationCondition.optimal:
        Optimal = True
    if Condition == pyo.TerminationCondition.maxTimeLimit or Condition == pyo.TerminationCondition.maxIterations:
        LimitStop = True
    if Optimal or LimitStop:
        try:
            WriteSolution = True
            Model.solutions.load_from(Results)   # Defer loading results until now, in case there is no solution to load
            SolverData = Results.Problem._list
            SolutionLB = SolverData[0].lower_bound
            SolutionUB = SolverData[0].upper_bound
        except:
            WriteSolution = False
            LimitStop = False
    return WriteSolution, LimitStop, SolutionLB, SolutionUB

Write output

Having processed the results to see what we've able to output, now we do the actual output, as shown in Figure 11.

The outputs are simple, consisting of some summary information about the model, the objective function value, and key results that are of interest – like the results shown in the Excel version of this model.

Figure 11. Write output

def Output(WriteSolution, LimitStop, SolutionLB, SolutionUB):
    print('Model: ', Model.name)
    print('Status:', Results.solver.termination_condition)
    print('Solver:', pyo.value(Model.Engine), '\n')
    if LimitStop:
        print('Objective bounds')
        print('----------------')
        if SolutionLB is None:
            print('Lower:      None')
        else:
            print(f'Lower: {SolutionLB:9,.2f}')
        if SolutionUB is None:
            print('Upper:      None\n')
        else:
            print(f'Upper: {SolutionUB:9,.2f}\n')
    if WriteSolution:
        print(f'Total cost = ${Model.TotalCost():11,.2f}\n')
        pd.set_option('display.max_rows', None)
        pd.set_option('display.max_columns', None)
        pd.options.display.float_format = '{:.0f}'.format

        Summary = pd.DataFrame()
        for p in Model.People:
            PeopleUsed = sum(pyo.value(Model.Allocation[s]) for s in Model.Shifts if pyo.value(Model.ShiftLabel[s]) == pyo.value(Model.PeopleLabel[p]))
            Summary.loc[pyo.value(Model.PeopleLabel[p]), 'Used'] = PeopleUsed
            Summary.loc[pyo.value(Model.PeopleLabel[p]), 'Available'] = pyo.value(Model.Available[p])
            Summary.loc[pyo.value(Model.PeopleLabel[p]), 'At least'] = pyo.value(Model.PeopleMin[p])
        display(Summary)

        Detail = pd.DataFrame()
        for s in Model.Shifts:
            if pyo.value(Model.Allocation[s]) >= 0.5:   # Use 0.5 to allow for precision issues
                for t in Model.Slots: 
                    Detail.loc[str(s) + '-' + pyo.value(Model.ShiftLabel[s]), str(pyo.value(Model.SlotLabel[t]))[0:5]] = pyo.value(Model.ShiftSlots[s, t]) * pyo.value(Model.Allocation[s])
        for s in Model.Shifts:
            for t in Model.Slots: 
                Detail.loc['Total', str(pyo.value(Model.SlotLabel[t]))[0:5]] = sum(pyo.value(Model.ShiftSlots[s, t] * Model.Allocation[s]) for s in Model.Shifts)
                Detail.loc['Required', str(pyo.value(Model.SlotLabel[t]))[0:5]] = pyo.value(Model.SlotRequired[t])
                Detail.loc['Surplus', str(pyo.value(Model.SlotLabel[t]))[0:5]] = sum(pyo.value(Model.ShiftSlots[s, t] * Model.Allocation[s]) for s in Model.Shifts) - pyo.value(Model.SlotRequired[t])
        display(Detail)
    else:
        print('No solution loaded')
        print('Model:')
        Model.pprint()

Main procedure

Finally, we bring everything together by calling each of the functions, starting with getting the data, and finishing with writing the output, as shown in Figure 12. Because the model is encapsulated within the global Model object, we don't need to pass it to the functions.

Figure 12. Main procedure

DefineData()
CheckData(Display=True)
DefineModel()
Results = SolveModel(Verbose=False)
WriteSolution, LimitStop, SolutionLB, SolutionUB = ProcessResults()
Output(WriteSolution, LimitStop, SolutionLB, SolutionUB)

Comparison between Excel and Python models

Optimal solution

The optimal solutions found by the Excel OpenSolver and Python Pyomo models are identical. Since this model has alternative optima, in general there's no guarantee that the two versions will find the same solution – even though both models use the CBC solver. That is, subtleties in how the model is presented to the solver may lead to it taking a different path during the solve process, resulting in a different alternative optima.

Model build process

The Python model has about 270 lines of code, though some of those lines are not strictly necessary (like def CheckData). Even so, the code appears complex – which it is, because we must manage every detail of the modelling process.

Conversely, the Excel model appears much simpler. That's largely because the spreadsheet takes care of some modelling aspects for us. But that appearance is somewhat deceptive. The Excel workbook contains more than 1,000 formulae (mostly in the Intermediate section, where we calculate the shift patterns multiplied by the usage variables). There are also quite complex relationships between the formulae, most of which is hidden from view behind the nice, neat grid interface.

Overall, the Excel model version was easier to build – mainly because it involves less low-level data manipulation for building the model and writing outputs. That's a key reason why Excel can be a useful tool for building prototype optimization models. An Excel model might even be sufficient for operational use. Though the flexibility of Python is generally better for models that require different data inputs. That's especially the case for models where the data varies in size, as that can be difficult to handle within the rigid grid structure of a spreadsheet.

Conclusion

In this article we replicated in Python Pyomo the model previously built in Excel OpenSolver. The two models produce the same result, but they are very different in terms of build process. Each model version has its advantages and disadvantages, as is typical when looking at the choice of tools for a modelling project.

If you would like to know more about this model, or you want help with your own models, then please contact us.

Schedule staff with enumerated shifts, Pyomo