本文介紹如何用數(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āng)庫(kù), 列代表客戶. 矩陣中的每一個(gè)值代表對(duì)應(yīng)的倉(cāng)庫(kù)到客戶的單位運(yùn)輸成本. 我們的目標(biāo)最小化總的運(yùn)輸成本.
下面我們用數(shù)學(xué)語言描述該問題.
輸入
- 倉(cāng)庫(kù)的供給量
,
, 其中
是倉(cāng)庫(kù)總數(shù)
- 客戶的需求量
,
, 其中
是客戶總數(shù)
- 倉(cāng)庫(kù)
到客戶
的單位運(yùn)輸成本是
輸出
- 需要計(jì)算倉(cāng)庫(kù)
到客戶
的運(yùn)輸量
下面我們寫出問題的目標(biāo)和約束.
目標(biāo)是最小化總的運(yùn)輸成本, 即
我們需要滿足的約束條件有兩個(gè):
- 每個(gè)倉(cāng)庫(kù)的出庫(kù)量不能超過其供給量:
,
- 每個(gè)客戶的需求應(yīng)該被滿足:
,
綜上所述, 我們可以把運(yùn)輸問題用線性規(guī)劃(Linear Programming)來表示.
標(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)有和
, 其中
代表倉(cāng)庫(kù),
代表客戶.
參數(shù)(Parameters)
參數(shù)是問題的輸入. 以上述運(yùn)輸問題為例, 我們的參數(shù)是: 供給量(), 需求量(
, 單位運(yùn)輸成本
.
決策變量(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()
數(shù)獨(dú)(Sudoku)
把數(shù)字1-9填入下圖的空格子中, 且滿足如下三個(gè)條件:
- 每個(gè)區(qū)塊 (圖中灰色方框包含的3
3小格子)包含數(shù)字1-9
- 每行包含數(shù)字1-9
- 每列包含數(shù)字1-9
我們通過數(shù)學(xué)規(guī)劃的方式求解該問題.
指標(biāo)
-
-- 填入的數(shù)字,
-
-- 第
行區(qū)塊, 區(qū)塊一共三行, 因此
-
-- 第
行區(qū)塊, 區(qū)塊一共三列, 因此
-
-- 區(qū)塊中元素的行, 每個(gè)區(qū)塊包含三行, 因此
-
-- 區(qū)塊中元素的列, 每個(gè)區(qū)塊包含三列
參數(shù)
-
-- 考慮第
行
列的區(qū)塊, 它的
行
列是否數(shù)字
決策變量
-
-- 考慮第
行
列的區(qū)塊, 它的
行
列是填入否數(shù)字
約束
- 已經(jīng)存在的值不能修改.
,
- 一個(gè)單元格同時(shí)只允許填入一個(gè)數(shù)字.
- 每個(gè)區(qū)塊包含數(shù)字1-9.
- 每行包含數(shù)字1-9.
- 每列包含數(shù)字1-9.
綜上所述, 我們的規(guī)劃可以寫成下面的整數(shù)規(guī)劃(Integer Programming). 注意: 無優(yōu)化目標(biāo).
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()