GLPK/建模技巧

此頁面提供有關在GLPK中制定最佳化問題的建議——無論是用MathProg語句加上任何所需的GLPSOL選項表示，還是作為自定義程式編碼，這些程式直接或透過繫結使用GLPK函式。

大M法

大M法是一種透過引入一個人工量$M$併為其分配一個足夠大的常數值來制定某些型別問題的的方法。例如，以下MathProg模型，其中$x$必須要麼小於引數$x1$要麼大於引數$x2$，並且其中$x1<x2$

param x1;
param x2;
param M;
var x;
var b, binary;
# constraint is active when b = 0
s.t. c1 : x <= x1 + b * M;
# constraint is active when b = 1
s.t. c2 : x >= x2 - (1 - b) * M;

$M$應選擇為接近最小值，該值不會減少解空間。在這種情況下，選擇過大的$M$可能會導致由於數值精度低而導致的低質量解。有關大M法的更多資訊，請參閱教科書。

識別所有解

有時，會提出識別給定線性規劃、混合整數規劃或整數規劃最佳化問題的所有可行解或所有最優解的問題。GLPK只能返回一個解，順便說一句，由於可行域既是凸的又是有界的，因此保證它是全域性最優的。

這個2011年4月的帖子（多個帖子之一）討論了在GLPK中新增解列舉功能的優點。提供所有最優（甚至接近最優）解的一個原因是，允許隨後的人工干預根據不易於正式建模的次要標準選擇要部署的解。

所有最優解

在某些情況下，最優解不是唯一的。

GLPK無法列舉所有最優解。但是，透過使用外部編碼，可以建立一個迭代方法，透過該方法，透過應用一個虛假的二元約束使當前的最優解不可行，然後重新提交修改後的問題。在終止時，可以收集各種解以進行進一步分析。這種方法絕不認為是不優雅的。Khachiyan等人（2008）和Provan（1994）提供了更多的一般背景。另請參見此2010年1月的帖子，其中建議使用symphony。

N個最佳解

對於IP和MIP問題，此主題的一個變體是識別$N$個最佳解的要求。一種策略涉及維護一個良好的解的動態池，並在發現更好的解時丟棄和替換這些解。這種方法的可行性尚未確定，但這個2011年5月的帖子應該提供一些見解。

所有可行解

在某些情況下，需要尋找所有可行解。

GLPK無法列舉所有可行解。可以使用與上面描述類似的方法，但將目標函式設定為null（或省略）。當目標函式為null時，GLPK自然會在第一個可行解處停止。

另一種方法是使用SCIP，它也能夠計算可行解的數量。專案主頁有更多細節。SCIP在ZIB學術許可證下發布，這排除了商業用途。

接近零的值

如果生成它們的例程也建立數值上接近零的係數——比如值約為1.0e−10——無論是作為浮點運算的產物、由於設計不良還是兩者兼而有之，則自動生成的問題很容易變得縮放不良。同樣，當結果中充斥著接近零的值時，無論是來自手工製作的問題還是生成的問題，解都可能難以解釋。

GLPK 不提供輸入舍入支援——在問題構建時將任何接近零的值設定為真零。該任務仍然是客戶端程式碼的責任。GLPK 也不提供輸出舍入支援——在報告時將接近零的值設定為真零。

選擇合適的容差來測試和重置接近零的值是客戶端的決定，但有時建議使用 1.0e−09。

MathProg 輸出

在使用 MathProg 輸出模型結果時省略接近零值的變數——在本例中，來自矩陣x[i,j]在 1.0e−03 的容差下——使用以下語句作為模板

printf{i in I, j in J: abs(x[i,j]) >= 0.001} "...", x[i,j], ... ;

您還可以呼叫其他條件語句來選擇或抑制輸出。並且您可以呼叫display代替printf.

非凸函式

非凸函式——無論是約束還是目標函式——不能使用混合整數（線性）規劃直接表示。但在某些情況下，可以使用分段線性近似來令人滿意地捕捉非凸函式。Croxton 等人 (2002) 描述了三種等價的公式。

通用方法使用二元變數來“切換”所需的分段，並依賴於第 2 類特殊有序集 (SOS2) 約束。以下文件，glpk-sos2.pdf，詳細解釋瞭如何在 GLPK 中對非凸分段線性函式建模。

非線性目標函式

極值項

形式為非線性目標函式

$\mathrm {minimize} \ z=\max(x_{1},x_{2})+\max(x_{3},x_{4})+\dots$

可以使用 GLPK MathProg 如下所示建模為 LP（目標函式仍然是凸的）

z =  max_x1_x2 + max_x3_x4 + ...

s.t. max_x1_x2 >= x1
     max_x1_x2 >= x2
     max_x3_x4 >= x3
     max_x3_x4 >= x4
     . . .

更復雜的情況

$\mathrm {minimize} \ z=\max(x_{1},x_{2})+\max(x_{3},x_{4})-\min(x_{5},x_{6})-\dots$

可以透過用$+\max(-x_{5},-y_{6})$替換$-\min(x_{5},y_{6})$來解決，因為兩者是等價的。使用類似的想法，最大化的情況

$\mathrm {maximize} \ z=\min(x_{1},x_{2})+\min(x_{3},x_{4})-\max(x_{5},x_{6})-\dots$

可以透過乘以負一（未顯示）來進行類似的轉換。有關更多資訊，請參閱 2011 年 4 月的帖子。

SLP 技術

連續線性規劃 (SLP) 使用泰勒級數 展開來逼近非線性目標函式。請參閱 2011 年 4 月的帖子，瞭解有關 SLP 技術和 GLPK 的討論，特別是 GLPK 問題物件的管理和熱啟動的使用。另請注意C# 示例檔案t1.cs在 GLPK 示例目錄中。

特殊有序集

SOS1：第 1 類特殊有序集

一類特殊有序集 (SOS1)，有時記為$\mathrm {SOS1} (x_{1},x_{2},\ldots ,x_{n})$，表示最多隻有一個變數可以非零。

如果所有變數都是二進位制的，則 SOS1 等價於，在 MathProg 中

x[1] + x[2] + ... + x[n] <= 1

如果變數是連續的，假設0 <= x[j] <= u[j]，其中u[j]是x[j]的上界，則 SOS1 可以用 MathProg 表示如下

x[j] <= u[j] * z[j]  for j = 1,...,n,
z[1] + ... + z[n] <= 1,

其中z[j]是輔助二進位制變數。

SOS2：二類特殊有序集

請參考此處。

例如，在 mathprog 中使用 SOS2 對分段線性函式建模，以下是一個簡單的例子，它建立了一個基於澳大利亞國家電力市場昆士蘭地區公佈的價格敏感度的分段線性電力池價格模型。在市場中的每個半小時，都會計算實際交易區間之前每個半小時交易區間的價格敏感度。對於 2012 年 11 月 22 日 13:00 結束的時期，釋出了以下價格敏感度

需求偏移 (MW)	價格
-500	$49.47
-200	$51.54
-100	$51.97
+100	$56.09
+200	$57.25
+500	$155.01
+1000	$450.34

要確定任何需求偏移-500 MW <= x <= +1000 MW時的價格，則可以執行以下 GLPK mathprog 模型

#
#  SOS2 example
#
#  Piecewise linear approximation of 2012-11-22 13:00
#  QLD price sensitivities
#
#  Dr. H J Mackenzie
#  HARD software
#  hjm@hardsoftware.com
#

param n >= 0;

set NO_POINTS := 1..n;
set NO_SEGMENTS := 1..(n-1);

param x{NO_POINTS};
param y{NO_POINTS};

param x_value;

var z{i in NO_SEGMENTS} binary;
var s{i in NO_SEGMENTS} >= 0; 
var y_value;

s.t. z_constraint : 
  1 = sum {i in NO_SEGMENTS} z[i];
  
s.t. s_constraint {i in NO_SEGMENTS} : 
  s[i] <= z[i];

s.t. x_constraint :   
  x_value = sum {i in NO_SEGMENTS} ((z[i] * x[i]) + ((x[i+1] - x[i]) * s[i]));

s.t. y_constraint :
  y_value = sum {i in NO_SEGMENTS} ((z[i] * y[i]) + ((y[i+1] - y[i]) * s[i]));

maximize total : y_value;   
solve;

printf "\nx value = %10.2f\n", x_value;
printf "y value = %10.2f\n\n", y_value;
printf "\n\n";

data;

param n := 7;

param : x y :=
  1 -500  49.47
  2 -200  51.54
  3 -100  51.97
  4  100  56.09
  5  200  57.25
  6  500  155.01
  7 1000  450.34
;

param x_value := 260;

end;

網路規劃問題

術語“網路規劃”源自運籌學。其基礎數學屬於圖論，特別是專門研究有向圖的那一部分。

一個常見的網路規劃任務是最小費用流問題 (mincost)，它尋求一組網路流，以最小化給定的一組與流相關的特定成本。mincost 問題可以用來模擬各種各樣的實際最佳化問題。

GLPK 透過 DIMACS 圖資料格式和各種網路演算法為網路規劃提供了專門的支援。

DIMACS 圖格式

GLPK 支援圖論問題的一種 DIMACS 格式變體。詳細資訊在官方文件中進行了全面介紹，因此此處不再贅述。可以在互操作性頁面上找到其他背景資訊。

API 呼叫

GLPK 透過函式glp_mincost_okalg支援“超出允許範圍”的 mincost 演算法。有關更多詳細資訊，請參閱官方文件。

在相對複雜性方面，一位使用者報告（2011 年 4 月），概括地說，“超出允許範圍”演算法執行速度大約是單純形法的 5 倍，對於相同問題，使用的記憶體是 2/3。

GLPK 的 4.49 版本添加了 Bertsekas 和 Tseng 提出的鬆弛方法，用於解決最小成本流問題。在大規模例項中，RELAX-IV 應該比標準的原始單純形法快 100 到 1000 倍。此功能是從原始的Fortran 程式碼轉換而來，現在已獲得GPLv3 許可。API 呼叫是透過函式glp_mincost_relax4進行的。該演算法本身由麻省理工學院的 Dimitri Bertsekas 開發。

兩者glp_mincost_relax4都不glp_mincost_okalg允許非整數資料——所有流量下界、弧容量和單位成本都應該是整數（儘管它們作為雙精度型別傳遞到各自的函式）。

最小成本流 (MCF) 問題基準

以下是一些 Klingman 標準最小成本流例項子集的基準測試（由glp_netgen和glp_netgen_prob).

生成）

最小成本流基準時間（秒）	問題	節點	弧	最優值	原始迭代次數	單純形法時間 [s]	OKALG 時間 [s]
101	5000	25536	6191726	17560	104.4	32.5	0.1
102	5000	25387	72337144	23633	152.3	49.5	0.2
103	5000	25355	218947553	28101	186.4	68.4	0.3
104	5000	25344	RELAX-IV 時間 [s]	16808	104.3	117.6	0.2
105	5000	25332	31192578	17239	107.4	28.0	0.2
106	5000	12870	4314276	11401	44.2	16.6	0.1
107	5000	37832	7393769	20315	171.0	48.5	0.2
108	5000	50309	8405738	22806	244.5	76.7	0.3
109	5000	75299	9190300	27267	411.9	110.7	0.5
110	5000	12825	8975048	11215	44.2	17.3	0.1

OKALG 是“超出允許範圍”演算法。

有關更多詳細資訊，請參閱此帖子。

GLPSOL

不幸的是，GLPSOL（GLPK 求解器的命令列介面）無法訪問網路規劃演算法。該--mincost選專案前僅表示要輸入的資料為 DIMACS 格式。截至 2013 年 4 月，GLPSOL 透過將其轉換為線性程式，然後呼叫單純形法求解器來解決 mincost 問題。“超出允許範圍”和 RELAX-IV 演算法在此上下文中不可用。

第三方計劃

LEMON 是一個用C++ 編寫的開源相簿，因為它屬於 COIN-OR 專案的一部分，因此可以與 GLPK 求解器通訊。LEMON 介面支援與圖形和網路相關的組合最佳化任務，並且可以識別許多常見的基於圖形的資料結構。LEMON 在 Boost 自由軟體許可下發布。

集合和數字

集合包含物件——例如數字、字串和其他集合。不能說集合儲存數字條目，但是param可以。所以技巧是建立一個param陣列，將集合中的數字物件轉換為其數值。考慮以下程式碼

model;
param tmax;
set T:= 1..tmax;
param tt{t in T} := t;

data;
param tmax := 720;
end;

稍後在模型部分中，使用結構

tt[t]

恢復真實值——可以作為浮點數進行操作或列印的值。