Benders分解由Jacques F. Benders在1962年提出^[1]. 它是一種把線性規(guī)劃問題分解為小規(guī)模子問題的技巧. 通過迭代求解主問題和子問題, 從而逼近原問題的最優(yōu)解. 與列生成(Column Generation)相比, Benders分解是一種行生成(Row Generation)技巧: 主問題的約束來自子問題的解. 本文介紹的內(nèi)容基于 Georrion和Graves^[2].

問題描述

考慮如下的數(shù)學(xué)規(guī)劃問題:

其中, , , . 和可以是非線性的. 可以是離散或連續(xù)的. 給定, 上述問題則變成了一個標(biāo)準(zhǔn)的線性規(guī)劃問題, 記作.

假設(shè). , 線性規(guī)劃存在最優(yōu)解. (否則原問題也不存在最優(yōu)解.）

Benders分解

等價形式

先把原問題寫成如下形式:

根據(jù)假設(shè), 問題

存在最優(yōu)解. 其對偶問題

的最優(yōu)解也存在. (不會寫對偶? 這篇文章手把手教你: <線性規(guī)劃技巧: 如何寫對偶問題>) .

因此進(jìn)而可以寫成:

寫成上述形式的好處是中的約束與變量無關(guān). 因此, 其最優(yōu)解一定是多面體(Polytope) 中的一個頂點. 令代表多面體頂點的集合. 因此原問題等價于

再把寫成數(shù)學(xué)規(guī)劃的形式:

問題分解

對原問題, 我們先固定得到, 然后考慮其對偶問題. 通過這種方式我們把原問題轉(zhuǎn)化成等價問題. 在新問題中, 每個約束條件對應(yīng)一個對偶問題的基本解, 且目標(biāo)函數(shù)與無關(guān). 這樣一來, 使得我們可以通過迭代的方式進(jìn)行求解: 主問題一開始只考慮少量的約束, 然后通過求解對偶問題添加新的約束到, 直到逼近最優(yōu)解.

主問題

我們知道是多面體的頂點. 是否需要把所有頂點加入到主問題中? 回顧主問題的目標(biāo)函數(shù). 顯然, 我們應(yīng)該找到使得最大化.

子問題

求解過程

思路

1. 初始化主問題. 令. 換句話說, 主問題一開始不考慮關(guān)于的約束. 令, 然后求解, 從而得到子問題的輸入.
2. 令OPT(master)代表最優(yōu)解的值. 由于主問題只考慮了部分約束, 因此OPT(master)是原問題最優(yōu)解的下界.
3. 求解子問題得到, 然后把約束添加到主問題.
4. 令OPT(sub)代表子問題最優(yōu)解的值. 回顧原問題的目標(biāo)函數(shù), 因此+OPT(sub)是原問題最優(yōu)解的上界.
5. 反復(fù)求解主問題和子問題直到上下界相等(或非常接近).

偽代碼

Init: B = empty; UB = inf; e = -1e-4;
solve master problem M(y, z) by fixing z := 0 and get y;
let LB := f(y) + z;
While UB - LB > e: 
    solve sub problem S(u|y) and get u;
    update UB := min {UB, f(y) + OPT(sub)};
    add constraint [b-F(y)]^T u <= z to the master problem;
    solve master problem M(y, z) and get y;
    update LB := OPT(master);
EndWhile

說明

1. 迭代過程中LB不會減小, UB不會增大.
2. 如果子問題的最優(yōu)解滿足唯一性, 這個迭代過程收斂(證明略).

例子: 設(shè)施選址問題(A Facility Location Problem)

某個工廠需要開設(shè)一些倉庫, 用來向它的客戶提供倉儲和配送服務(wù).

• 考慮個候選倉庫和個客戶.
• 新建一個倉庫有開倉成本. 服務(wù)一個客戶有連接成本, 例如為客戶提供配送服務(wù)和退換貨服務(wù)帶來的成本等.
• 我們需要決定開設(shè)哪些倉庫, 并指定每個客戶對應(yīng)的服務(wù)倉庫.
• 目標(biāo)是最小化開倉成本與所有客戶的連接成本之和.

圖片來自網(wǎng)絡(luò).

(如上圖所示, 藍(lán)色方框代表開設(shè)的倉庫, 其它候選倉庫未展示. 紅色圓圈代表客戶. 黑色的線代表服務(wù)關(guān)系.)

整數(shù)線性規(guī)劃

Indices

• -- 倉庫
• -- 客戶

Parameters

• -- 倉庫的開倉成本
• -- 倉庫服務(wù)客戶的連接成本

Decision Variables

• -- 倉庫是否服務(wù)客戶
• -- 倉庫是否開設(shè)

如果開設(shè)倉庫確定, 我們指定離客戶最近的倉庫來服務(wù)這個客戶, 從而達(dá)到連接成本最低. 因此, 設(shè)施選址問題的核心是決定開設(shè)哪些倉庫, 因而上述規(guī)劃的決策變量的取值范圍可以寫成. 我們得到如下規(guī)劃:

Benders分解

給定, 的最優(yōu)解等價于如下問題的最優(yōu)解.

它的對偶問題則是我們要求解的子問題:

結(jié)合原問題的目標(biāo), 我們得到主問題:

注意: 我們需要檢查本文開頭提到的假設(shè)是否成立. 即, 對任意可行的, 子問題是否存在最優(yōu)解. 事實上, 當(dāng)時, 子問題的目標(biāo)函數(shù)的最大值是, 因而不存在最優(yōu)解! 因此我們需要為主問題添加約束來保證假設(shè)條件成立. 注意到原問題的最優(yōu)解至少要保證開設(shè)1個倉庫, 把它加入主問題的約束從而保證了子問題最優(yōu)解的存在(想想為什么).

主問題

說明 主問題的約束條件一開始為空, 即. 它的約束條件通過求解子問題得到.

Python實現(xiàn)

主問題模型

# master_model.py

from ortools.linear_solver import pywraplp

import numpy as np


class MasterModel(object):

    def __init__(self, f, fix_z=None):
        """
        :param f: 開倉成本(m維向量)
        :param fix_z: 限定z的值. 目的是初始化的時候令z=0,然后求解主問題得到y(tǒng)的值.
        """
        self._solver = pywraplp.Solver('MasterModel',
                                       pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
        self._f = f
        self._fix_z = fix_z
        self._m = len(self._f)
        self._y = None  # 決策變量
        self._z = None  # 決策變量
        self._solution_y = None  # 計算結(jié)果
        self._solution_z = None  # 計算結(jié)果
        # 迭代求解之前會調(diào)用add_constraint
        # 所以決策變量的初始化要放在前面
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()

    def _init_decision_variables(self):
        self._y = [self._solver.NumVar(0, 1, 'y[%d]' % i)
                   for i in range(self._m)]
        self._z = self._solver.NumVar(-self._solver.Infinity(), self._solver.Infinity(), 'z')

    def _init_objective(self):
        obj = self._solver.Objective()
        for i in range(self._m):
            obj.SetCoefficient(self._y[i], self._f[i])
        obj.SetCoefficient(self._z, 1)
        obj.SetMinimization()

    def add_constraint(self, alpha, beta):
        """ Benders主流程會調(diào)用此方法為主問題增加新的約束.

        :param alpha: 子問題的解
        :param beta: 子問題的解
        """
        ct = self._solver.Constraint(sum(alpha), self._solver.Infinity())
        for i in range(self._m):
            ct.SetCoefficient(self._y[i], sum(beta[i]))
        ct.SetCoefficient(self._z, 1)

    def _init_constraints(self):
        if self._fix_z is not None:
            ct = self._solver.Constraint(self._fix_z, self._fix_z)
            ct.SetCoefficient(self._z, 1)

        ct = self._solver.Constraint(1, self._solver.infinity())
        for i in range(self._m):
            ct.SetCoefficient(self._y[i], 1)

    def solve(self):
        self._solver.Solve()
        self._solution_y = [self._y[i].solution_value() for i in range(self._m)]
        self._solution_z = self._z.solution_value()

    def get_solution(self):
        return self._solution_y, self._solution_z

    def get_objective_value(self):
        return sum(np.array(self._solution_y) * np.array(self._f)) + self._solution_z

子問題模型

# sub_model.py

from ortools.linear_solver import pywraplp
import numpy as np


class SubModel(object):

    def __init__(self, y, C):
        """
        :param y: 主問題的解(代表y_i=1代表開設(shè)倉i)
        :param C: 連接成本(m*n維矩陣)
        """
        self._solver = pywraplp.Solver('SubModel',
                                       pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
        self._y = y
        self._c = C
        self._m = len(self._c)
        self._n = len(self._c[0])
        self._alpha = None  # 決策變量
        self._beta = None  # 決策變量
        self._constraints = []  # 所有約束(后續(xù)獲取對偶變量, 得到原始問題的解)
        self._solution_alpha = None  # 計算結(jié)果
        self._solution_beta = None  # 計算結(jié)果

    def _init_decision_variables(self):
        self._alpha = [self._solver.NumVar(-self._solver.Infinity(),
                                           self._solver.Infinity(),
                                           'alpha[%d]' % j)
                       for j in range(self._n)]

        self._beta = [[self._solver.NumVar(0, self._solver.Infinity(), 'beta[%d][%d]' % (i, j))
                      for j in range(self._n)] for i in range(self._m)]

    def _init_constraints(self):
        for i in range(self._m):
            self._constraints.append([])
            for j in range(self._n):
                ct = self._solver.Constraint(-self._solver.Infinity(), self._c[i][j])
                ct.SetCoefficient(self._alpha[j], 1)
                ct.SetCoefficient(self._beta[i][j], -1)
                self._constraints[i].append(ct)

    def _init_objective(self):
        obj = self._solver.Objective()
        for j in range(self._n):
            obj.SetCoefficient(self._alpha[j], 1)
        for i in range(self._m):
            for j in range(self._n):
                obj.SetCoefficient(self._beta[i][j], -self._y[i])

        obj.SetMaximization()

    def solve(self):
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()
        self._solver.Solve()
        self._solution_alpha = [self._alpha[j].solution_value() for j in range(self._n)]
        self._solution_beta = [[self._beta[i][j].solution_value()
                                for j in range(self._n)]
                               for i in range(self._m)]

    def get_solution(self):
        return self._solution_alpha, self._solution_beta

    def get_objective_value(self):
        return sum(self._solution_alpha) - sum(np.array(self._y) * np.sum(self._solution_beta, axis=1))

    def get_dual_values(self):
        """ 得到原問題的解: x[i][j]
        """
        return [[self._constraints[i][j].dual_value()
                for j in range(self._n)]
                for i in range(self._m)]

Benders分解流程

# benders_proc.py

import numpy as np

from master_model import MasterModel
from sub_model import SubModel


class BendersProc(object):
    """ Benders分解流程
    """

    def __init__(self, f, C, max_iter=10000):
        """
        :param f: 開倉成本(m維向量)
        :param C: 連接成本(m*n矩陣), 其中m是候選設(shè)施的數(shù)量, n是客戶的數(shù)量
        :param max_iter: 最大循環(huán)次數(shù)
        """
        self._f = f
        self._c = C
        self._iter_times = 0
        self._max_iter = max_iter
        self._status = -1  # -1:執(zhí)行錯誤; 0:最優(yōu)解; 1: 達(dá)到最大循環(huán)次數(shù)
        self._ub = np.inf
        self._lb = -np.inf
        self._solution_x = None  # 計算結(jié)果
        self._solution_y = None  # 計算結(jié)果

    def _stop_criteria_is_satisfied(self):
        """ 根據(jù)上下界判斷是否停止迭代
        """
        if self._ub - self._lb < 0.0001:
            self._status = 0
            return True
        if self._iter_times >= self._max_iter:
            if self._status == -1:
                self._status = 1
            return True
        return False

    def solve(self):
        # 初始令z=0. 求解主問題得到y(tǒng)
        master0 = MasterModel(self._f, fix_z=0)
        master0.solve()
        y, z = master0.get_solution()
        # 下面的迭代需要重新生成master對象
        # 因為master0中z=0是約束條件
        master = MasterModel(self._f)
        sub = None
        # 迭代過程
        while not self._stop_criteria_is_satisfied():
            # 求解子問題
            sub = SubModel(y, self._c)
            sub.solve()
            # 更新上界
            fy = sum(np.array(self._f) * np.array(y))
            self._ub = min(self._ub, fy + sub.get_objective_value())
            # 生成主問題的約束
            alpha, beta = sub.get_solution()
            master.add_constraint(alpha, beta)
            master.solve()
            # 更新y和下界
            y, z = master.get_solution()
            self._lb = master.get_objective_value()

            print(">>> iter %d: lb = %.4f, ub = %.4f" % (self._iter_times, self._lb, self._ub))
            self._iter_times += 1

        # 保存結(jié)果
        self._solution_x = sub.get_dual_values()
        self._solution_y = y

        status_str = {-1: "error", 0: "optimal", 1: "attain max iteration"}
        print(">>> Terminated. Status:", status_str[self._status])

    def print_info(self):
        print("---- Solution ----")
        res = {}
        for i in range(len(self._f)):
            if self._solution_y[i] == 0:
                continue
            connected_cities = [str(j) for j in range(len(self._c[0]))
                                if self._solution_x[i][j] > 0]
            res[i] = ', '.join(connected_cities)

        for f, c in res.items():
            print("open facility: %d, connected cities: %s" % (f, c))

主函數(shù)

# main.py

from benders_proc import BendersProc
from data import f, C  # data.py包含了一個實例


if __name__ == '__main__':
    bp = BendersProc(f, C)
    bp.solve()
    bp.print_info()

完整代碼

參考文獻(xiàn)

1. J.F. Benders. Partitioning Procedures for Solving Mixed-Variables Programming Problems. Numerische Mathematik, Vol. 4, 1962. ?

2. A.M. Geoffrion and G.W. Graves. Multicommodity distribution system design by benders decomposition. Management Science 20, no. 5, 822–844, 22, 1974. ?