Fix and Optimize mit abstrakten Unfixierern

Blaupause für ein leicht adaptierbares Fix-and-Optimize Framework.

Ausgangslage

In diesem Artikel entwerfe ich ein simples Fix-and-Optimize Framework mit Hilfe der Google OR-Tools und abstrakten Fixier-Klassen.

Fix-and-Optimize

Egal welches Problem mit Fix-and-Optimize gelöst werden soll, der Ansatz ist immer sehr ähnlich:

Die Variablen (i.d.R. die binären- und/oder ganzzahligen Variablen) eines mathematischen Optimierungsmodells werden in zwei Teilmengen aufgeteilt:

• Eine Menge von Variablen deren Werte fixiert sind (z.B. auf 0 oder 1) und die vom Solver nicht verändert werden dürfen
• Eine kleine Menge von Variablen, deren Wert nicht fixiert sind und deren Werte (unter berücksichtigung der fixierten Variablen) vom Solver bestimmt werden sollen.

Dies führt dazu, dass immer nur eine kleine Anzahl von Entscheidungen durch den Solver getroffen werden müssen und das korrespondierende mathematische Optimierungsmodell meist sehr schnell gelöst werden kann. Nachteil dieses Vorgehens ist, dass die Lösungen nicht zwingend optimal sind.

erzeuge Startlösung
while Abbruchkriterium nicht erfüllt:

    löse die Fixierung einiger Variablen (gemäß eines Musters)
    löse das resultierende MIP

    if Lösung ist besser als beste Bekannte:
        aktualisiere beste Lösung
    
    fixiere alle Variablen auf ggf. neue Werte

Wir müssen also sicherstellen, dass drei zentrale Konzepte implementiert werden:

• Ein mathematisches Optimierungsmodell
• Eine Fixierungsroutine
• Eine Routine, die gemäß einer Logik Variablenwerte löst
• Ein iteratives Vorgehen, das fixiert und löst und Lösungen updated.

Die konkreten Ausprägungen der beiden Konzepte sind stark problemspezifisch. So fixieren wir in eine Standortplanungsproblem eher eröffnete Standorte, während in einem Losgrößenplanungsproblem Rüstperioden fixiert werden.

(Abstraktes) Mathematisches Optimierungsmodell

Jede Fix-and-Optimize Heuristik benötigt ein mathematisches Optimierungsmodell. Unabhängig davon, welches Modell gewählt wird, sind für das Verfahren zwei Funktionen relevant:

• das Lösen
• das Abrufen des Zielfunktionswerts

Deswegen wird die abstrakte Modellklasse wie folgt entworfen:

class MipModel(ABC):
    def __init__(self,name="mip") -> None:
        self.solver = pywraplp.Solver(name,pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
    
    def solve(self):
        self.solver.Solve()
    
    def objective(self):
        return self.solver.Objective().Value()

Beispielimplementierung: MIP SLULSP

class Slulsp(MipModel):

    def __init__(self,d, s, l) -> None:

        self.solver = pywraplp.Solver("SLULSP",pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
        self.periods = len(d)

        # Variablen
        self.L = self.__initL(self.solver, self.periods)
        self.x = self.__initX(self.solver, self.periods)
        self.y = self.__initY(self.solver, self.periods)

        # Nebenbedingungen
        self.__storageConst(self.solver,self.L, self.x,self.periods,d)
        self.__setupConst(self.solver,self.x, self.y, self.periods,d)

        # ZF
        self.__addObjective(self.solver,l,self.L,s, self.y)
    
    
    def __initL(self, solver, T):
        # Lagerbestand am Ende von t
        return [solver.NumVar(0.0, solver.infinity(), "") for _ in range(T)]   
    def __initX(self, solver, T):
        # Produktionsmenge X in t
        return [solver.NumVar(0.0, solver.infinity(), "") for _ in range(T)]   
    def __initY(self, solver, T):
        # = 1 wenn in t produziert wird, 0 sonst
        return [solver.BoolVar("") for _ in range(T)]                          

    def __storageConst(self,solver,L,x,T,d):
        solver.Add(x[0] - d[0] == L[0])
        for t in range(1,T):
            solver.Add(L[t-1] + x[t] - d[t] == L[t])

    def __setupConst(self,solver,x,y,T, d):
        # Hinreichend große Zahl
        M = sum(d)                      
        for t in range(T):
            solver.Add(x[t] <= M*y[t])

    def __addObjective(self,solver,l,L,s,y ):
        solver.Minimize( solver.Sum([l*L[t] + s*y[t] for t in range(T)]) )

Fixierung und unfixierung von Variablen: Generelle Idee

Wenn die Heuristik von vornherein einigermaßen effizient implementiert werden soll, dann betrachten wir in jeder Iteration nur all jene Variablen deren Fixierung wir lösen wollen. Alle anderen Variablen des Modells lassen wir unangetastet. e

Wir fixieren- bzw. lösen diese Variablen indem wir festlegen, dass ihre Lowerbound der Upperbound entspricht. Wollen wir beispielsweise festlegen, dass eine Binärvariable auf 0 fixiert werden soll, dann nutzen wir variable.SetBounds(0,0), um die Lowerbound (erste 0) und Upperbounds (zweite 0) auf Null zu setzen. Soll die Variable auf 1 fixiert werden, schreiben wir variable.SetBounds(1,1).

Dieses Vorgehen erledigt was wir wollen, es gilt aber etwas wichtiges zu beachten:

Sobald die Lower- bzw. Upperbound irgendeiner Variable im Modell verändert wird, werden alle aktuellen Lösungswerte, die die Variablen haben zurückgesetzt.

Dieser Umstand (den ich schmerzlich lernen musste, Grüße gehen raus an den Debugger), wird dann problematisch, wenn wir nacheinander die Bounds einer Variable auf den Wert einer aktuellen Lösung fixieren wollen. Denn mit er Fixierung der ersten Variable, sind die weiteren Werte nicht mehr abrufbar.

Beispiel anhand von SLULSP (crasht)

mip = Slulsp(d,s,l)
mip.solve()

for y in mip.y:
    # Crasht ab der 2. Iteration
    wert = y.solution_value()
    y.SetBounds(wert,wert)

Aus diesem Grund gehen wir wie folgt vor:

• zunächt speichern wir eine Referenz auf die entsprechende Variable sowie deren aktuellen Wert in einer Klasse zwischen
• nach dem vollständigen Zwischenspeichern setzen wir die Bounds

Hierfür nutzen wir die folgende Zwischenspeicherklasse:

@dataclass
class VarValueStorage:
    variable : object
    value : float

Nachfolgend nutzen wir das beigefügte Muster, damit das Programm sauber durchläuft.

Beispiel anhand von SLULSP (crasht nicht mehr)

mip = Slulsp(d,s,l)
mip.solve()

zwischenspeicher = []
# Speichert alle Werte
for y in mip.y:
    wert = y.solution_value()
    zwischenspeicher.append(VarValueStorage(y, wert))

# Setzt die Bounds ohne zu crashen
for zw_speicher in zwischenspeicher: 
    zw_speicher.variable.SetBounds(zw_speicher.value, zw_speicher.value)

Abstrakte Unfixer/Fixer Klasse

Mit dem zuvor dargestellten Ansatz entwerfen wir jetzt ein Unfixer Klasse. Hier werden je nach Problemstellung unterschiedliche Ansätze implementiert, um Variablen zu lösen.

class Unfixer(ABC):

    def __init__(self) -> None:
        self.variables_to_unfix = []
    
    @abstractmethod
    def unfix(self,model : MipModel):
        """
        Löst die Fixierung von Variablen. 
        Welche Variablen konkret ausgewählt werden, hängt von der Implementierung ab.
        Wichtig: self.variables_to_unfix leeren.
        """
        pass

    def update_fix_values(self):
        """ Fixiert die neuen (besseren) Werte der Variablen, die zuvor gelöst wurden.
        """

        # Aktualisiert die neuen Werte der Variablen
        for stored_var in self.variables_to_unfix:
            new_value = stored_var.variable.solution_value()
            stored_var.value = new_value
    
    def refix(self):
        """ Fixiert die Variablen auf die Werte, die sie vor der Unfixierung hatten """
        for stored_var in self.variables_to_unfix:
            prev_value = stored_var.value
            stored_var.variable.SetBounds(prev_value,prev_value)

Beispielimplementierung: Unfixer für SLULSP

Nachfolgend wollen wir rollierend immer lookahead aufeinanderfolgende Rüstperioden unfixieren. Wenn wir über den Planungshorizont hinaus lappen, beginnen wir vorne.

Wir lösen daher die Fixierung der y-Variablen, so dass sie wieder die Werte 0 oder 1 annehmen können. Um den aktuellen fixierten Wert zwischenzuspeichern reicht es aus, die aktuelle Lower Bound zu speichern (da sie ja der Upperbound entspricht).

Da SlulspUnfixer die abstrakte Klasse Unfixer implementiert, brauchen wir uns über update_fix_values und refix keine Sorgen machen. Diese sind direkt verfügbar.

class SlulspUnfixer(Unfixer):

    def __init__(self, lookahead : int) -> None:
        super().__init__()
        self.t_now = 0
        self.lookahead = lookahead

    def unfix(self, model: Slulsp):

        self.variables_to_unfix = []
        
        for pseudo_t in range(self.t_now,self.t_now+self.lookahead):
            t = pseudo_t % len(model.y)
            stor = VarValueStorage(model.y[t], model.y[t].lb())
            self.variables_to_unfix.append(stor)
        
        for var in self.variables_to_unfix:
            var.variable.SetBounds(0,1)
        
        self.t_now = self.t_now + 1 % len(model.y)

Fix-and-Optimize Framework

Zu guter letzt wird das Fix and Optimize Framework implementiert, das nach einer definierten Anzahl Iterationen abbricht.

class FixAndOptimize:

    def __init__(self,model : MipModel, unfixer, minimize:bool, iterations=100) -> None:
        self.model = model
        self.unfixer = unfixer
        self.iterations = iterations
        self.minimize = minimize
    
    def __is_better(self,val_a : float, val_b : float):
        if self.minimize:
            return val_a < val_b
        return val_a > val_b
    
    def solve(self):
        """ Setzt voraus, dass eine initiale Lösung gegeben und fixiert ist."""

        # Speichere Startzielfunktionswert als Besten
        self.bestObjective = self.model.objective()

        # Fix and Optimize für eine gegebene Anzahl Iterationen
        for i in range(self.iterations):

            # Löst Variablenwerte
            self.unfixer.unfix(self.model)

            # Optimiert
            self.model.solve()

            # Updatet die beste bekannte Lösung, falls zutreffend
            if self.__is_better(self.model.objective(), self.bestObjective):
                self.bestObjective = self.model.objective()
                self.unfixer.update_fix_values()
            
            # Fixiert die zuvor gelösten und reoptimierten Variablen erneut
            self.unfixer.refix()

Konkrete Implementierung: Fix-and-Optimize

# Nachfrage der Perioden
d = [50,30,20,10,5,50,10,30]    
# Rüstkostensatz
s = 50                          
# Lagerkostensatz
l = 1                           


# Startlösung bei der in jeder Periode gerüstet wird
mip = Slulsp(d,s,l)
for v in mip.y:
    v.SetBounds(1,1)
mip.solve()

# Unfixer mit zwei aufeinanderfolgenden Perioden
unfixer = SlulspUnfixer(lookahead=2)
f_and_o = FixAndOptimize(mip, unfixer,True)
f_and_o.solve()

Try yourself

Der oben dargestellte Code liegt auf GitHub.

Ausprobiert werden kann er auf Google Colab