參考:python運籌優化(一):Cplex for python使用簡介
下面是一個簡單的優化模型:
$$ min \sum_{i=1}^{n} \sum_{j=1}^{m} c_{ij} x_{ij} $$
$$s.t.$$
$$\sum_{i=1}^{n} a_{ij}x_{ij} \le b_j \qquad \forall j$$
$$x_{ij} \ge l_{ij} \qquad \forall i,j$$
$$ x_{ij} \le u_{ij} \qquad \forall i,j$$
在上述優化例子中,n、m、a、c、l、u、b為輸入參數,假設為給定。為了編寫Python代碼,我們將這些參數設置如下:
有10*5個決策變量$x_{ij}$,10行5列
# 初始化變量的值 import random n = 10 m = 5 set_I = range(1, n+1) set_J = range(1, m+1)
后面的部分
# 初始化已知量 # 通過 random 生成隨機數 # 都是以 字典 的形式來生成的 # 按照公式中需要的數據一一生成 c = {(i, j): random.normalvariate(0, 1) for i in set_I for j in set_J} a = {(i, j): random.normalvariate(0, 5) for i in set_I for j in set_J} l = {(i, j): random.randint(0, 10) for i in set_I for j in set_J} u = {(i, j): random.randint(10, 20) for i in set_I for j in set_J} b = {j: random.randint(0, 30) for j in set_J}
1. 導入運籌優化庫
# 導入運籌優化庫 import docplex.mp.model as cpx opt_model = cpx.Model(name="MIP Model")
2. 定義決策變量
在這一步之后,我們有一個名為opt_model的模型對象。接下來,我們需要添加決策變量。在Python字典(或panda系列)中存儲決策變量是標准的,其中字典鍵是決策變量,值是決策變量對象。一個決策變量由三個主要屬性定義:它的類型(連續、二進制或整數)、它的下界(默認為0)和上界(默認為無窮大)。對於上面的例子,我們可以將決策變量定義為:
# 如果 x 是連續的話,按照下面的定義 x_vars_c = {(i, j): opt_model.continuous_var(lb=l[i, j], ub=u[i, j], name="x_{0}_{1}".format(i, j)) for i in set_I for j in set_J}
不同的定義方式,針對上面的問題,下面兩個只是參考
# 如果 x 是二進制 x_vars_b = {(i, j): opt_model.binary_var(name="x_{0}_{1}".format(i, j)) for i in set_I for j in set_J} # 如果 x 是整數 x_vars = {(i, j): opt_model.integer_var(lb=l[i, j], ub=u[i, j], name="x_{0}_{1}".format(i, j)) for i in set_I for j in set_J}
3. 約束條件
在設置決策變量並將它們添加到我們的模型之后,就到了設置約束的時候了。任何約束都有三個部分:左手邊(通常是決策變量的線性組合)、右手邊(通常是數值)和意義(小於或等於、等於、大於或等於)。要設置任何約束,我們需要設置每個部分:
$$\sum_{i=1}^{n} a_{ij}x_{ij} \le b_j \qquad \forall j$$
# <= constraints, 小於等於 # 也是寫成了字典的形式 # 這樣可以方便后面查看調用 # 實際上使用 opt_model.add_constraint 就可以將其加入到模型內部了 constraints = {j: opt_model.add_constraint( ct=opt_model.sum(a[i,j] * x_vars_c[i,j] for i in set_I) <= b[j], ctname="constraint_{0}".format(j)) for j in set_J}
下面兩個作為參考
# >= constraints, 大約等於 constraints = {j: opt_model.add_constraint( ct=opt_model.sum(a[i,j] * x_vars_c[i,j] for i in set_I) >= b[j], ctname="constraint_{0}".format(j)) for j in set_J} # == constraints, 等於 constraints = {j: opt_model.add_constraint( ct=opt_model.sum(a[i,j] * x_vars_c[i,j] for i in set_I) == b[j], ctname="constraint_{0}".format(j)) for j in set_J}
下面是另外兩個邊界的約束條件
$$x_{ij} \ge l_{ij} \qquad \forall i,j$$
constraints_l = {(i, j): opt_model.add_constraint( ct=x_vars_c[i,j] >= l[i,j], ctname="constraint_l_{0}_{1}".format(i,j)) for i in set_I for j in set_J}
$$ x_{ij} \le u_{ij} \qquad \forall i,j$$
constraints_u = {(i, j): opt_model.add_constraint( ct=x_vars_c[i,j] <= u[i,j], ctname="constraint_u_{0}_{1}".format(i,j)) for i in set_I for j in set_J}
4. 目標函數
下一步是定義一個目標,它是一個線性表達式。我們可以這樣定義目標:
$$ min \sum_{i=1}^{n} \sum_{j=1}^{m} c_{ij} x_{ij} $$
# 兩層循環 # 然后將所有的數據相加 objective = opt_model.sum(x_vars_c[i,j] * c[i,j] for i in set_I for j in set_J) # for maximization opt_model.maximize(objective)
作為參考的
# for minimization opt_model.minimize(objective)
5. 求解模型
# solving with local cplex opt_model.solve()
作為參考的
# solving with cplex cloud opt_model.solve(url="your_cplex_cloud_url", key="your_api_key")
6. 獲得結果
現在我們完成了。我們只需要得到結果並進行后期處理。panda包是一個很好的數據處理庫。如果問題得到最優解,我們可以得到和處理結果如下:
# 可以直接打印結果 opt_model.print_solution()
輸出結果
objective: 216.829 x_1_1=2.000 x_1_2=11.000 x_1_3=6.000 x_1_4=17.000 x_1_5=12.000 x_2_2=5.000 x_2_3=17.000 x_2_4=12.000 x_2_5=8.000 x_3_1=12.000 x_3_2=5.000 x_3_3=14.000 x_3_4=10.000 x_3_5=8.000 x_4_1=2.000 x_4_3=4.000 x_4_4=14.000 x_4_5=2.000 x_5_2=16.321 x_5_3=5.000 x_5_4=1.000 x_5_5=19.000 x_6_1=12.000 x_6_2=9.000 x_6_3=5.000 x_6_4=11.000 x_6_5=10.000 x_7_1=5.000 x_7_2=3.000 x_7_3=1.000 x_7_4=19.000 x_7_5=13.895 x_8_1=13.000 x_8_2=17.000 x_8_3=10.000 x_8_4=10.000 x_8_5=4.000 x_9_1=18.000 x_9_2=2.000 x_9_3=5.000 x_9_4=10.000 x_9_5=16.000 x_10_1=19.000 x_10_2=4.000 x_10_3=18.000 x_10_4=10.000 x_10_5=5.000
通過 pandas.DataFrame 顯示
# 合並在一起的效果 import pandas as pd # 將 x_var_c 變量加入到 dataframe 里面 opt_df = pd.DataFrame.from_dict(x_vars_c, orient="index", columns = ["variable_object"]) opt_index = pd.MultiIndex.from_tuples(opt_df.index, names=["colums_i", "colums_j"]) # 生成了新的索引 opt_df.reset_index(inplace=True) # CPLEX # variable_object 顯示的是 x_1_1,實際上是一個 object,因此可以調用響應的函數 # 對於整個列進行操作 # 獲取每個對象的最優結果 opt_df["solution_value"] = opt_df["variable_object"].apply(lambda item: item.solution_value) opt_df.drop(columns=["variable_object"], inplace=True) opt_df.to_csv("./optimazation_solution.csv")
這里,opt_df是一個包含每個決策變量$x_{ij}$的最優值的panda dataframe。我們還可以將這些結果保存到CSV文件中,如上所示。
我們只討論了Python中的高級建模,但是上面的所有包都包含有用的函數和數據結構,在編寫准備生產的代碼時應該考慮這些函數和數據結構。例如,在gu中,可以使用opt_model.addVars()一次性添加一組變量,而在CPLEX中是opt_model.continuous_var_dict()、opt_model.binary_var_dict()或opt_model.integer_var_dict(),在PuLP中可以使用plp.LpVariable.dicts()。
7. 代碼合並在一起,去掉多余的部分
# 初始化變量的值 import random n = 10 m = 5 set_I = range(1, n+1) set_J = range(1, m+1) # 通過 random 生成隨機數,都是以 字典 的形式來生成的,按照公式中需要的數據一一生成 c = {(i, j): random.normalvariate(0, 1) for i in set_I for j in set_J} a = {(i, j): random.normalvariate(0, 5) for i in set_I for j in set_J} l = {(i, j): random.randint(0, 10) for i in set_I for j in set_J} u = {(i, j): random.randint(10, 20) for i in set_I for j in set_J} b = {j: random.randint(0, 30) for j in set_J} # 導入運籌優化庫 import docplex.mp.model as cpx opt_model = cpx.Model(name="MIP Model") # 如果 x 是連續的話,按照下面的定義 x_vars = {(i, j): opt_model.continuous_var(lb=l[i, j], ub=u[i, j], name="x_{0}_{1}".format(i, j)) for i in set_I for j in set_J} # 3個constraints,約束條件 for j in set_J: opt_model.add_constraint(ct=opt_model.sum(a[i,j] * x_vars[i,j] for i in set_I) <= b[j]) for i in set_I: for j in set_J: opt_model.add_constraint(ct=x_vars[i,j] >= l[i,j]) opt_model.add_constraint(ct=x_vars[i,j] <= u[i,j]) # 兩層循環,然后將所有的數據相加 objective = opt_model.sum(x_vars[i,j] * c[i,j] for i in set_I for j in set_J) # for maximization opt_model.maximize(objective) # solving with local cplex opt_model.solve() # 可以直接打印結果 opt_model.print_solution()