用數(shù)學(xué)規(guī)劃的方式求解優(yōu)化問題

本文介紹如何用數(shù)學(xué)語言對(duì)實(shí)際中的優(yōu)化問題進(jìn)行建模. 通過建立數(shù)學(xué)模型, 我們利用現(xiàn)成的求解器可以便捷地計(jì)算出最優(yōu)解(或可行解).

運(yùn)輸問題

考慮三個(gè)糧食儲(chǔ)量分別是100, 200, 300的倉(cāng)庫(kù) (單位:噸, 下文省略). 我們需要把糧食運(yùn)送給4個(gè)客戶, 其需求分別是: 120, 60, 270, 150.

倉(cāng)庫(kù)到客戶的單位運(yùn)輸成本用矩陣C描述:
\begin{align} \begin{bmatrix} 350 & 200 & 300 & 250 \\ 220 & 330 & 300 & 270 \\ 215 & 230 & 290 & 240 \\ \end{bmatrix} \end{align}
其中行代表倉(cāng)庫(kù), 列代表客戶. 矩陣中的每一個(gè)值代表對(duì)應(yīng)的倉(cāng)庫(kù)到客戶的單位運(yùn)輸成本. 我們的目標(biāo)最小化總的運(yùn)輸成本.

下面我們用數(shù)學(xué)語言描述該問題.

輸入

  • 倉(cāng)庫(kù)的供給量s_i, i=1, 2, ... ,m, 其中m是倉(cāng)庫(kù)總數(shù)
  • 客戶的需求量d_j, j= 1, 2, ..., n, 其中n是客戶總數(shù)
  • 倉(cāng)庫(kù)i到客戶j的單位運(yùn)輸成本是c_{i, j}

輸出

  • 需要計(jì)算倉(cāng)庫(kù)i到客戶j的運(yùn)輸量x_{i,j}

下面我們寫出問題的目標(biāo)和約束.

目標(biāo)是最小化總的運(yùn)輸成本, 即
\min \sum_{i,j}c_{ij}x_{ij}.

我們需要滿足的約束條件有兩個(gè):

  • 每個(gè)倉(cāng)庫(kù)的出庫(kù)量不能超過其供給量: \sum_{j} x_{ij} \leq a_i, \forall i
  • 每個(gè)客戶的需求應(yīng)該被滿足: \sum_{i} x_{ij} = d_j, \forall j

綜上所述, 我們可以把運(yùn)輸問題用線性規(guī)劃(Linear Programming)來表示.

\begin{aligned} \min~& \sum_{ij}c_{ij} x_{ij} \\ \text{s.t. } & \sum_{j} x_{ij} \leq a_i, \forall i \\ & \sum_{i} x_{ij} = d_j, \forall j \\ & x_{ij} \geq 0, \forall i, j. \end{aligned}

標(biāo)準(zhǔn)實(shí)踐

為了更加直觀地寫出數(shù)學(xué)模型, 我們可以總結(jié)一份標(biāo)準(zhǔn)的指南. 它包含四個(gè)基本步驟:

  • 指標(biāo)(Indices)
    指標(biāo)的作用是主要為了簡(jiǎn)化記號(hào). 以上述運(yùn)輸問題為例, 我們的指標(biāo)有ij, 其中i代表倉(cāng)庫(kù), j代表客戶.

  • 參數(shù)(Parameters)
    參數(shù)是問題的輸入. 以上述運(yùn)輸問題為例, 我們的參數(shù)是: 供給量(s_i), 需求量(d_j), 單位運(yùn)輸成本c_{i,j}.

  • 決策變量(Decision Variables)
    決策變量是算法的輸出.

  • 優(yōu)化目標(biāo)(Objective)
    一般是最小化或最大化一個(gè)目標(biāo)函數(shù). 在某些情況下, 問題只需要找到一個(gè)可行解, 因此也可以不指定優(yōu)化目標(biāo).

  • 約束(Constraints)
    用等式或不等式描述解的限制.

求解規(guī)劃

常用的商用求解器有Gruobi和CPLEX(可申請(qǐng)教育和學(xué)術(shù)的lisense). 商用求解器功能強(qiáng)大, 能求解多種類型的規(guī)劃問題, 例如整數(shù)規(guī)劃, 混合整數(shù)規(guī)劃, 二次規(guī)劃等. 免費(fèi)的求解器有Google的ORtools, 它把一些開源的求解器做了集成, 求解速度雖然比不上商用求解器, 實(shí)際中也能滿足很多業(yè)務(wù)需求.

求解方式有兩種:

第一種是直接用商用求解器提供的IDE. 按照求解器的建模語法把模型寫出來, 然后求解. 建模語法的好處是非常貼近公式化的描述, 所見即所得.

第二種是調(diào)用求解器提供的API, 初始化參數(shù), 約束, 目標(biāo), 然后求解.

本文我們使用開源工具ORtools求解(基本的教程請(qǐng)自行g(shù)oogle,需要翻墻)

Python實(shí)現(xiàn)

模型

# model.py

from ortools.linear_solver import pywraplp
import numpy as np


class TransportModel(object):

    def __init__(self, a, d, C):
        """
        :param a: 供給量(m維向量), m代表倉(cāng)庫(kù)數(shù)量
        :param d: 需求量(n維向量), n代表客戶數(shù)量
        :param C: 單位運(yùn)輸成本(m*n維矩陣), C[i][j]代表倉(cāng)庫(kù)i到客戶j的單位運(yùn)輸成本
        """
        self._solver = pywraplp.Solver('TransportModel',
                                       pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
        self._a = a
        self._d = d
        self._C = C
        self._m = len(self._a)  # 倉(cāng)庫(kù)數(shù)量
        self._n = len(self._d)  # 客戶數(shù)量
        self._x = None  # 決策變量
        self._solution_x = None  # 計(jì)算結(jié)果
        self._obj_val = None  # 目標(biāo)函數(shù)值

    def _init_decision_variables(self):
        self._x = [
            # 0 <= x[i][j] <= infinity
            [self._solver.NumVar(0, self._solver.infinity(), "x[%d][%d]" % (i, j))
             for j in range(self._n)] for i in range(self._m)
        ]

    def _init_constraints(self):
        # 每個(gè)倉(cāng)庫(kù)的出庫(kù)量不能超過其供給量
        # sum(x[i][j]) <= a[i], over j
        for i in range(self._m):
            ct = self._solver.Constraint(0, self._a[i])
            for j in range(self._n):
                ct.SetCoefficient(self._x[i][j], 1)
        # 每個(gè)客戶的需求應(yīng)該被滿足
        # sum(x[i][j]) == b[j], over i
        for j in range(self._n):
            ct = self._solver.Constraint(self._d[j], self._d[j])
            for i in range(self._m):
                ct.SetCoefficient(self._x[i][j], 1)

    def _init_objective(self):
        obj = self._solver.Objective()
        for i in range(self._m):
            for j in range(self._n):
                obj.SetCoefficient(self._x[i][j], self._C[i][j])
        obj.SetMinimization()

    def solve(self):
        self._init_decision_variables()
        self._init_constraints()
        self._init_objective()
        self._solver.Solve()
        # 求解器返回的解
        self._solution_x = [
            [self._x[i][j].solution_value() for j in range(self._n)]
            for i in range(self._m)
        ]
        # sum(C[i][i] * x[i][j]) over i,j
        self._obj_val = np.sum(np.array(self._C) * np.array(self._solution_x))

    def print_result(self):
        print("最優(yōu)值 = ", self._obj_val)
        print("最優(yōu)解 x = ")
        print(np.array(self._solution_x))

主函數(shù)

# main.py

from data import a, d, C  # 運(yùn)輸問題實(shí)例
from model import TransportModel


if __name__ == '__main__':
    tm = TransportModel(a, d, C)
    tm.solve()
    tm.print_result()

完整代碼: 運(yùn)輸問題

數(shù)獨(dú)(Sudoku)

把數(shù)字1-9填入下圖的空格子中, 且滿足如下三個(gè)條件:

  1. 每個(gè)區(qū)塊 (圖中灰色方框包含的3\times3小格子)包含數(shù)字1-9
  2. 每行包含數(shù)字1-9
  3. 每列包含數(shù)字1-9

我們通過數(shù)學(xué)規(guī)劃的方式求解該問題.

指標(biāo)

  • n -- 填入的數(shù)字, n \in \{1, 2, ..., 9 \}
  • i -- 第i行區(qū)塊, 區(qū)塊一共三行, 因此i \in \{1, 2, 3\}
  • j -- 第j行區(qū)塊, 區(qū)塊一共三列, 因此j \in \{1, 2, 3\}
  • p -- 區(qū)塊中元素的行, 每個(gè)區(qū)塊包含三行, 因此p \in \{1,2,3 \}
  • q -- 區(qū)塊中元素的列, 每個(gè)區(qū)塊包含三列q \in \{ 1, 2, 3\}

參數(shù)

  • a_{i,j,p,q,n} \in \{ 0, 1\} -- 考慮第ij列的區(qū)塊, 它的ij列是否數(shù)字n

決策變量

  • x_{i,j,p,q,n} \in \{ 0, 1\} -- 考慮第ij列的區(qū)塊, 它的ij列是填入否數(shù)字n

約束

  1. 已經(jīng)存在的值不能修改.
    x_{i,j,p,q, n} \geq a_{i,j,p, q, n}, \forall i,j,p,q, n
  2. 一個(gè)單元格同時(shí)只允許填入一個(gè)數(shù)字.
    \sum_n x_{i,j,p,q,n} = 1, \forall i,j,p,q
  3. 每個(gè)區(qū)塊包含數(shù)字1-9.
    \sum_{p, q} x_{i,j, p, q, n} = 1, \forall i, j, n
  4. 每行包含數(shù)字1-9.
    \sum_{j, q} x_{i,j,p,q,n} = 1, \forall i,p, n
  5. 每列包含數(shù)字1-9.
    \sum_{i, p} x_{i,j,p,q,n} = 1, \forall j,q, n

綜上所述, 我們的規(guī)劃可以寫成下面的整數(shù)規(guī)劃(Integer Programming). 注意: 無優(yōu)化目標(biāo).

\begin{aligned} \min~& 0 \\ \text{s.t. } & x_{i,j,p,q,n} \geq a_{i,j,p,q,n}, \forall i, j,p,q, n \\ & \sum_n x_{i,j,p,q,n} = 1, \forall i,j,p,q \\ & \sum_{p, q} x_{i,j, p, q, n} = 1, \forall i, j, n \\ & \sum_{j, q} x_{i,j,p,q,n} = 1, \forall i,p, n \\ & \sum_{i, p} x_{i,j,p,q,n} = 1, \forall j,q, n \\ & x_{i,j,p,q} \in \{ 0,1\} . \end{aligned}

Python實(shí)現(xiàn)

模型

# model.py

from ortools.linear_solver import pywraplp
import numpy as np


class SudokuModel(object):

    def __init__(self, a):
        """
        :param a: Sudoku實(shí)例 
        """
        self._solver = pywraplp.Solver('SudokuModel',
                                       pywraplp.Solver.BOP_INTEGER_PROGRAMMING)
        self._a = a
        self._x = None  # 決策變量
        self._solution_x = None  # 計(jì)算結(jié)果

    def __init_decision_variables(self):
        self._x = np.empty((3, 3, 3, 3, 9)).tolist()
        for i in range(3):
            for j in range(3):
                for p in range(3):
                    for q in range(3):
                        for n in range(9):
                            # 已知數(shù)字不允許修改
                            # x[i][j][p][q][n] >= a[i][j][p][q][n]
                            self._x[i][j][p][q][n] \
                                = self._solver.IntVar(self._a[i][j][p][q][n], 1,
                                                      'x[%d][%d][%d][%d][%d]' % (i, j, p, q, n))

    def __init_constraints(self):
        # 一個(gè)單元格同時(shí)只允許填入一個(gè)數(shù)字
        # sum(x[i][j][p][q][n]) = 1, over n
        for i in range(3):
            for j in range(3):
                for p in range(3):
                    for q in range(3):
                        ct = self._solver.Constraint(1, 1)
                        for n in range(9):
                            ct.SetCoefficient(self._x[i][j][p][q][n], 1)
        # 每個(gè)區(qū)塊包含數(shù)字1-9
        # sum(x[i][j][p][q][n]) = 1, over p, q
        for i in range(3):
            for j in range(3):
                for n in range(9):
                    ct = self._solver.Constraint(1, 1)
                    for p in range(3):
                        for q in range(3):
                            ct.SetCoefficient(self._x[i][j][p][q][n], 1)
        # 每行包含數(shù)字1-9
        # sum(x[i][j][p][q][n]) = 1, over j, q
        for i in range(3):
            for p in range(3):
                for n in range(9):
                    ct = self._solver.Constraint(1, 1)
                    for j in range(3):
                        for q in range(3):
                            ct.SetCoefficient(self._x[i][j][p][q][n], 1)
        # 每列包含數(shù)字1-9
        # sum(x[i][j][p][q][n]) = 1, over i, p
        for j in range(3):
            for q in range(3):
                for n in range(9):
                    ct = self._solver.Constraint(1, 1)
                    for i in range(3):
                        for p in range(3):
                            ct.SetCoefficient(self._x[i][j][p][q][n], 1)

    def solve(self):
        self.__init_decision_variables()
        self.__init_constraints()
        self._solver.Solve()
        self._get_solution_x()

    def _get_solution_x(self):
        self._solution_x = np.empty((3, 3, 3, 3))
        for i in range(3):
            for j in range(3):
                for p in range(3):
                    for q in range(3):
                        for n in range(9):
                            if self._x[i][j][p][q][n].solution_value() == 1:
                                self._solution_x[i][j][p][q] = n + 1

    def print_result(self):
        res = np.empty((9, 9))
        for i in range(3):
            for p in range(3):
                for j in range(3):
                    for q in range(3):
                        res[i*3+p][j*3+q] = self._solution_x[i][j][p][q]
        print(res)

主函數(shù)

# main.py

from model import SudokuModel
from data import a


if __name__ == '__main__':
    sm = SudokuModel(a)
    sm.solve()
    sm.print_result()

完整代碼: Sudoku

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡(jiǎn)書系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 今天我們將用幾個(gè)綜合實(shí)例來實(shí)踐Lingo。 例1 求解非線性方程組 規(guī)劃問題本來就是給出優(yōu)化條件和限制條件,之后得...
    熱愛學(xué)習(xí)的高老板閱讀 1,702評(píng)論 0 3
  • 久違的晴天,家長(zhǎng)會(huì)。 家長(zhǎng)大會(huì)開好到教室時(shí),離放學(xué)已經(jīng)沒多少時(shí)間了。班主任說已經(jīng)安排了三個(gè)家長(zhǎng)分享經(jīng)驗(yàn)。 放學(xué)鈴聲...
    飄雪兒5閱讀 7,865評(píng)論 16 22
  • 今天感恩節(jié)哎,感謝一直在我身邊的親朋好友。感恩相遇!感恩不離不棄。 中午開了第一次的黨會(huì),身份的轉(zhuǎn)變要...
    余生動(dòng)聽閱讀 10,911評(píng)論 0 11
  • 可愛進(jìn)取,孤獨(dú)成精。努力飛翔,天堂翱翔。戰(zhàn)爭(zhēng)美好,孤獨(dú)進(jìn)取。膽大飛翔,成就輝煌。努力進(jìn)取,遙望,和諧家園??蓯塾巫?..
    趙原野閱讀 3,539評(píng)論 1 1
  • 在妖界我有個(gè)名頭叫胡百曉,無論是何事,只要找到胡百曉即可有解決的辦法。因?yàn)槭侵缓偞蠹乙杂瀭饔灲形摇皟A城百曉”,...
    貓九0110閱讀 3,728評(píng)論 7 3

友情鏈接更多精彩內(nèi)容