PMCM-2019-F題-問題1（3）
- 作業要求
- 進一步剪枝
- 建立數學模型
- 繪制優化結果圖
- 建模及優化結果
- 以后寫博文會在同名公眾號“ManTou饅頭”上同步更新啦~ 主要原因是：
（1）博文容易被搬運
（2）博客回復大家的問題不及時(當然大家博客有問題還是可以回復和私信一起討論！ManTou會盡量一起討論解決的！)
（3）CSDN廣告多，影響觀感，公眾號沒有粉絲就不用擔心有廣告植
入了哦~- 代碼要結合前幾次作業的內容，需要整合一下~
- 每周的事情比較多，數學模型作業的質量可能不能讓大家滿意~ 大家多理解支持吖~
- 公眾號回復“PMCM2019F”，獲取源程式
拜托不要直接copy，拜托了拜托了
作業要求
- 進一步剪枝
在 作業（3）的基礎上，進一步減少網路鄰接矩陣中的有向邊的數量， - 建立模型
以 航跡長度盡可能小 和 經過矯正區域進行校正的次數盡可能少 為目標建立模型，
進一步剪枝
- 第二次剪枝（ PMCM-2019-F題-問題1（2）中的剪枝稱為第一次剪枝）
- 剪枝的標準
不走回頭路：洗掉方向和前進方向相反的有向邊， - 采用的方法
向量乘積小于零則說明是鈍角，也就說明是回頭路，這樣的邊就剪掉 - 代碼實作
# 第二次剪枝
AB = [x_ab[1] - x_ab[0], y_ab[1] - y_ab[0], z_ab[1] - z_ab[0]]
edge1 = copy.copy(edge) # 建立一個新的edge
for (i, j) in edge1:
if A_B[0] * (x_coord[j] - x_coord[i]) \
+ A_B[1] * (y_coord[j] - y_coord[i]) \
+ A_B[1] * (z_coord[j] - z_coord[i]) <= 0:
edge.remove((i, j))
print("第二次剪枝之后的邊數：", len(edge))
- 第三次剪枝
- 剪枝的標準
轉彎飛行的波動不太大：刪去距離AB路線太遠的校正點，以AB為軸，R為半徑的圓柱之外的點 - 采用的方法
計算校正點到直線的距離（參考鏈接），大于R的剪去， - 代碼實作
# 第三次剪枝
# 剪去以AB為軸，R為半徑的圓柱之外的點
# https://zhuanlan.zhihu.com/p/24760577
R = 10000
third_cut = []
for i in range(1, shape - 1):
if dist[i, shape - 1] * np.sqrt(1 - (((x_coord[shape - 1] - x_coord[i]) * (x_coord[shape - 1] - x_coord[0])
+ (y_coord[shape - 1] - y_coord[i]) * (y_coord[shape - 1] - y_coord[0])
+ (z_coord[shape - 1] - z_coord[i]) * (z_coord[shape - 1] - z_coord[0]))
/ (dist[0, shape - 1] * dist[i, shape - 1])) ** 2) >= R:
third_cut.append(i)
edge2 = copy.copy(edge) # 建立一個新的edge
for (i, j) in edge2:
if i in third_cut or j in third_cut:
edge.remove((i, j))
print("第三次剪枝之后的邊數：", len(edge))
- 剪枝的結果
剪枝可以有效減少計算量，可能對最終結果沒有太大的影響，
建立數學模型
- 模型的基本認識
- 模型應該是最短路0-1整數規劃模型，
x
i
j
{x_{ij}}
xij?是一個0-1變數，A點的入度為0出度為1，B點的出度為0入度為1，其他點入度等于出度（有向圖的出度和入度）,K是頂點項的懲罰系數（懲罰函式的作用）
min ? ∑ ( i , j ) ∈ E w i j x i j + K ∑ ( i , j ) ∈ E x i j \min \sum_{(i, j) \in E} w_{i j} x_{i j}+K \sum_{(i, j) \in E} x_{i j} min∑(i,j)∈E?wij?xij?+K∑(i,j)∈E?xij?
s.t. { ∑ ( i , j ) ∈ E x i j = ∑ ( j , k ) ∈ E x j k , ? j ≠ A , B ∑ ( A , j ) ∈ E x A j = 1 , ∑ ( j , A ) ∈ E x j A = 0 ∑ ( j , B ) ∈ E x j B = 1 , ∑ ( B , j ) ∈ E x B j = 0 x i j = 0 or 1 , ? i , j \left\{\begin{array}{l}\sum_{(i, j) \in E} x_{i j}=\sum_{(j, k) \in E} x_{j k}, \forall j \neq A, B \\ \sum_{(A, j) \in E} x_{A j}=1, \quad \sum_{(j, A) \in E} x_{j A}=0 \\ \sum_{(j, B) \in E} x_{j B}=1, \sum_{(B, j) \in E} x_{B j}=0 \\ x_{i j}=0 \text { or } 1, \quad \forall i, j\end{array}\right. ????????∑(i,j)∈E?xij?=∑(j,k)∈E?xjk?,?j?=A,B∑(A,j)∈E?xAj?=1,∑(j,A)∈E?xjA?=0∑(j,B)∈E?xjB?=1,∑(B,j)∈E?xBj?=0xij?=0 or 1,?i,j? - 水平和垂直校正點誤差不超過上限
v i ≤ α 1 , h i ≤ α 2 , i ∈ V v_{i} \leq \alpha_{1}, h_{i} \leq \alpha_{2}, i \in V vi?≤α1?,hi?≤α2?,i∈V
v i ≤ β 1 , h i ≤ β 2 , i ∈ H v_{i} \leq \beta_{1}, h_{i} \leq \beta_{2}, i \in H vi?≤β1?,hi?≤β2?,i∈H
v B ≤ θ , h B ≤ θ v_{B} \leq \theta, h_{B} \leq \theta vB?≤θ,hB?≤θ - 相鄰校正點誤差之間的關系
( I i h i + δ w i j ? h j ) x i j = 0 , ( i , j ) ∈ E \left(I_{i} h_{i}+\delta w_{i j}-h_{j}\right) \quad x_{i j}=0,(i, j) \in E (Ii?hi?+δwij??hj?)xij?=0,(i,j)∈E
( ( 1 ? I i ) v i + δ w i j ? v j ) x i j = 0 , ( i , j ) ∈ E \left(\left(1-I_{i}\right) v_{i}+\delta w_{i j}-v_{j}\right) x_{i j}=0,(i, j) \in E ((1?Ii?)vi?+δwij??vj?)xij?=0,(i,j)∈E - 采用大M法將限制條件線性化（大M法）
I i h i + δ w i j ? h j ≤ ( 1 ? x i j ) M , ( i , j ) ∈ E I_{i} h_{i}+\delta w_{i j}-h_{j} \leq\left(1-x_{i j}\right) M,(i, j) \in E Ii?hi?+δwij??hj?≤(1?xij?)M,(i,j)∈E
( 1 ? I i ) v i + δ w i j ? v j ≤ ( 1 ? x i j ) M , ( i , j ) ∈ E \left(1-I_{i}\right) v_{i}+\delta w_{i j}-v_{j} \leq\left(1-x_{i j}\right) M,(i, j) \in E (1?Ii?)vi?+δwij??vj?≤(1?xij?)M,(i,j)∈E
- Gurobi：Pyhon API
具體函式的使用方法，參考我的這篇博文：簡單使用 Gurobi：Python API - 基本步驟
- 添加變數 x
- 添加變數h和v
- 添加關鍵的約束
- 添加目標函式
- 代碼實作
# 建立模型
model1 = Model()
# 添加新變數x[i,j], i=0 to shape-1 ,j = 0 to shape-1
x = model1.addVars(shape, shape, vtype=GRB.BINARY, name='x')
# 添加新變數h[i],v[i],i = 0 to shape-1
h = model1.addVars(shape, vtype=GRB.CONTINUOUS, name='h')
v = model1.addVars(shape, vtype=GRB.CONTINUOUS, name='v')
# 添加限制條件 x[i,j]==0 根據矩陣資訊進行剪枝
for i in range(shape):
for j in range(shape):
if C[i, j] == 0:
model1.addConstr(x[i, j] == 0)
# 添加限制條件的起始點，中間點，終點的出度和入度條件
# 對于有向圖來說，頂點的出邊條數稱為該頂點的出度，頂點的入邊條數稱為該頂點的入度
out_degree = [0] * shape
in_degree = [0] * shape
# out_degree[i]表示第i個點的出度
for (i, j) in edge:
out_degree[i] = out_degree[i] + x[i, j]
# in_degree[i]表示第i個點的入度
for (j, i) in edge:
in_degree[i] = in_degree[i] + x[j, i]
for i in range(shape):
if i == 0:
# 起點A的出度為1，入度為0
model1.addConstr(out_degree[i] == 1)
model1.addConstr(in_degree[i] == 0)
elif 0 < i < shape - 1:
# 中間點的出度等于入度
model1.addConstr(out_degree[i] == in_degree[i])
else:
# 終點B的入度為1，出度為0
model1.addConstr(out_degree[i] == 0)
model1.addConstr(in_degree[i] == 1)
# 添加限制條件h，v進行約束，M = 100000
M = 10000
for (i, j) in edge:
model1.addConstr(prop[i] * h[i] + delta * dist[i, j] - h[j] <= M * (1 - x[i, j]))
model1.addConstr((1 - prop[i]) * v[i] + delta * dist[i, j] - v[j] <= M * (1 - x[i, j]))
for i in range(shape):
if i == 0:
# 起點的垂直和水平誤差為0
model1.addConstr(h[i] == 0)
model1.addConstr(v[i] == 0)
elif 0 < i < shape - 1:
# 中間點的誤差約束條件
if prop[i] == 1:
# 垂直校正點前的誤差約束條件
model1.addConstr(h[i] <= alpha2)
model1.addConstr(v[i] <= alpha1)
else:
# 水平校正點前的誤差約束條件
model1.addConstr(h[i] <= beta2)
model1.addConstr(v[i] <= beta1)
else:
# 終點前的垂直和水平誤差約束條件
model1.addConstr(h[i] <= theta)
model1.addConstr(v[i] <= theta)
# 添加目標函式
K = 10000
dicct_modify = {}
for i in range(shape):
for j in range(shape):
dicct_modify[i, j] = K
dicct_modify = tupledict(dicct_modify)
model1.setObjective(dict_dist.prod(x) + dicct_modify.prod(x), GRB.MINIMIZE)
model1.update()
model1.optimize()
for s in model1.getVars():
if s.x == 1:
print('%s' '%g' % (s.varName, s.x))
print('Obj:%g' % model1.objVal)
繪制優化結果圖
# 將優化的路徑點取出繪圖
node = [0, 503, 294, 91, 607, 540, 250, 340, 277, 612]
# 定義空陣列用來存放校正點
x_coord_node = []
y_coord_node = []
z_coord_node = []
for j in node:
x_coord_node.append(x_coord[j])
y_coord_node.append(y_coord[j])
z_coord_node.append(z_coord[j])
ax.plot(x_coord_node, y_coord_node, z_coord_node, linewidth=2, c='r')
# 保存繪制的圖形
plt.savefig("FigData1.jpg", bbox_inches='tight')
plt.show()
建模及優化結果
結果如圖所示~
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/200944.html
標籤：其他
上一篇：Python-心電預處理