數學規划模型

數學規划模型是運籌學的一個分支，用來研究在一定約束條件下如何求解目標函數的最優值。目標函數在約束條件下的極值問題。(s.t = subject to 約束於)

數學規划模型常分為以下幾類（不互斥）：

1. 線性規划：目標函數和約束條件是決策變量的線性函數，常用單純形法求解
2. 非線性規划：目標函數或約束條件是決策變量的非線性函數，沒有通用的解法，搜索求解
3. 整數規划：決策變量只取整數的數學規划
4. 0-1 規划：整數規划的特例

下面分別介紹各個模型，並且拓展了兩個常用模型：最大最小化模型和多目標規划模型

• 線性規划的求解

  無論是用 Matlab 還是 Python 求解，都需要先將模型轉化為標准型：

  \[min\quad C^T*X \\[0.3cm] s.t = \begin{cases} AX\le b \\[0.2cm] AeqX = beq \\[0.2cm] lb \le X \le ub \end{cases} \]

  其中，\(X\) 表示決策變量，\(C\) 表示決策變量系數，\(A\) 表示不等式（\(\le\)）約束的系數矩陣，\(Aeq\) 表示等式約束的系數矩陣，\(lb, ub\) 分別表示決策變量的上下界。（是最小化目標函數；不等式約束是小於等於約束）

  如下線性規划模型：

  \[min\quad Z = 4x_1 + 15x_2 + x_3 + 125x_4 \\[0.3cm] s.t = \begin{cases} 3x_1 + 3x_2 + 15x_4 \geq 32 \\[0.2cm] 5x_1 + 2x_3 + x_4 = 24 \\[0.2cm] 14x_1 + 7x_4 \leq 42 \\[0.2cm] x_1,x_2,x_4 \geq 0 \end{cases} \]

  轉化為標准型后，各參數值為：

  \(X = \left[\begin{matrix}x_1\\x_2\\x_3\\x_4\end{matrix}\right]\quad C = \left[\begin{matrix}4\\15\\1\\125\end{matrix}\right] \quad A = \left[\begin{matrix} -3&-3&0&-15\\ 14&0&0&7\end{matrix}\right]\quad b = \left[\begin{matrix}-32\\ 42\end{matrix}\right]\\ Aeq =\left[\begin{matrix}5&0&2&4\end{matrix}\right]\quad beq =\left[\begin{matrix}24\end{matrix}\right]\quad lb=\left[\begin{matrix}0\\0\\-\infty\\0\end{matrix}\right]\quad ub=\left[\begin{matrix}\infty\\\infty\\\infty\\\infty\end{matrix}\right]\)

  然后輸入如下命令求解（Python）：

  from scipy import optimize
  import numpy as np
  
  # 若各參數若沒有，就為空[]
  ans = optimize.linprog(C, A, b, Aeq, beq, lb, ub, X0, OPTION) # X0 表示 初始值矩陣
  
  

  若是 Matlab，則為：

  % fval 表示函數值，其他意義同上
  [X, fval] = linprog[C, A, b, Aeq, beq, lb, ub, X0]
  

  注意：若原模型求的最大值，其結果應該加負號，無窮多解時，只返回一個

  如上面的 \(x_1\) 不再是 \(x_1\geq 0\)，而是 \(x_1 > 0\)，那么上面的下界應該怎么寫呢？可以寫成 \(x_1 \geq 0.001\)，或者其他的數

• 線性整數規划的求解

  線性整數規划與線性規划相比，只是多了一個整數約束

  % intlinprog 不能指定初始值，intcon 參數是行向量，表示哪些決策變量是整數，如 intcon = [1, 3]
  [x, fval] = intlinprog[C, intcon, A, b, Aeq, beq, 1b, ub]
  

• 0-1線性規划

  0-1線性規划是線性整數規划的特例，依然可以用 intlinprog​ 函數，只不過需要在上下界做做手腳

  令 \(lb=\left[\begin{matrix}0\\0\\-\infty\\0\end{matrix}\right]\quad ub=\left[\begin{matrix}1\\1\\\infty\\1\end{matrix}\right]\) 表示 \(x_1, x_2, x_4\) 只取整數 \(0\ or\ 1\)，\(x_3\) 不做限制

• 非線性規划的求解

  同線性規划一樣，要先將非線性規划模型轉化為標准型：

  \[\begin{align} &min\quad f(x) \\[0.3cm] &s.t = \begin{cases} AX\leq b\quad AeqX = beq &線性約束 \\[0.2cm] C(x)\leq 0\quad Ceq(x) = 0 &非線性約束 \\[0.2cm] lb\leq x\leq ub \end{cases} \end{align} \]

  如下非線性規划模型：

  \[\begin{align} &min f(x) = x_1^2 + x_2^2 + x_3^2 + 8 \\[0.3cm] &s.t = \begin{cases} x_1^2 -x_2 + x_3^2 \geq 0\\[0.2cm] x_1 +x_2^2+x_3^2\leq 20\\[0.2cm] -x_1 - x_2^2 + 2 = 0\\[0.2cm] x_2 + 2x_3^2 = 3\\[0.2cm] x_1, x_2, x_3 \geq 0 \end{cases} \end{align} \]

  轉化為非線性標志型后，各參數為：

  \(A\ =\ b\ =\ Aeq\ =\ beq\ =\ ub\ =\ \left[\begin{matrix}\ \end{matrix}\right]\quad lb = \left[\begin{matrix}0\\0\\0\end{matrix}\right]\\C =\left[\begin{matrix} -x_1^2 + x_2 -x_3^2;\ x_1 + x_2^2 + x_3^2 - 20\end{matrix}\right]\\[0.4cm] Ceq = \left[\begin{matrix}-x_1 -x_2^2 + 2;\ x_2 + 2x_3^2 - 3\end{matrix}\right]\)

  Matlab 求解代碼為：

  [x, fval] = fmincon(@fun, X0, A, b, Aeq, beq, lb, ub, @nonlfun, option)
  

  注意事項：

  1. 非線性規划中初始值是必須的，大多情況下求解的是一個局部最優解。

  2. 如果要使求得的解更接近全局最優，兩種方法：

     • 選取多組初始值，在求得的最后結果中選取最優
     • 先用蒙特卡洛模擬（詳見蒙特卡洛模型專題），得到一個或多個蒙特卡洛解，在用這個（些）解作為初始值

  3. option 選項用於給定解法：內點法（默認，interior-point）、序列二次規划（sqp）、有效集法（active-set）、信賴域反射算法（trust-legion-reflective），不同算法有其適用場景，可以都試一遍選取最優

  4. @fun 表示目標函數，需要單獨編寫一個 m 文件（名字為fun）存儲目標函數，針對該題，如下：

     function f = fun(x)
     	f = x(1)^2 + x(2)^2 + x(3)^2 + 8
     end
     

  5. @nonlfun 表示非線性部分的約束，同樣需要編寫一個 m 文件存儲（名字為nonlfun），針對該題，如下：

     function [C, Ceq] = nonlfun(x)
     	C = [
     		-x(1)^2 + x(2) - x(3)^2;
     		x(1) + x(2)^2 + x(3)^2 - 20;
     	]
     	Ceq = [
     		-x(1) - x(2)^2 + 2;
     		x(2) + 2*x(3)^2 - 3;
     	]
     end
     

• 拓展之最大最小化模型

  何為最大最小化模型？如下題：

  給定多個城鎮的坐標和一個區域，問在該區域內何處建立急救中心，使得急救中心離各城鎮的距離中的最大值最小？

  這是一個簡單的問題，則可以寫得目標函數為：$\min \limits_{(x,\ y)}\ {\max\limits_{i}\ [\ |x -a_i| + |y - b_i|\ ]} $，假設求得是曼哈頓距離

  [X, fval] = fminmax(@Fun, X0, A, b, Aeq, beq, lb, ub, @nonlfun, option)
  

  這里的 @Fun 也是一個 m 文件，寫法如下：

  function f = Fun(x)
  a = [];
  b = []; % 各城鎮的坐標
  f = zeros(m, 1); % 多個目標函數
  for i = 1: m
  	f(i) = abs(x(1) - a(i)) + abs(x(2) - b(i)); % 曼哈頓距離
  end
  end
  

• 多目標規划模型

  多目標規划常常通過加權組合轉化為單目標規划模型，有以下注意點：

  1. 需要將每個目標函數統一化為求最大值或最小值。（通過加負號）

  2. 加權組合時需要考慮到量綱（單位）問題，若量綱不同，則需要標准化去量綱操作，所謂的標准化有兩個作用：去掉單位（通過 \(\frac{單位}{單位}\)）；將數據都變成同一數量級。常用的標准化方法有：

     \[\begin{align} 最大-最小值規范化\quad \quad &X_{new} = \frac{X_{old}\ -\ X_{min}}{X_{max}\ -\ X_{min}}\\[0.3cm] 均值-方差規范化\quad \quad&X_{new} = \frac{X_{old} - \overline{X}}{S}\\[0.3cm] &S = \sqrt{\frac{1}{n-1}\sum_{i = 1}^{n}(X_i - \overline{X})^2}\\[0.3cm] 歸一化方法\quad \quad&X_{new} = \frac{X_{old}}{\sum_{i = 1}^{n}X_i} \end{align} \]

     當然，當題目給出了每個目標函數相應的參考值時，也可以通過 \(\frac{f}{參考值}\) 來標准化。

  3. 權重問題

  4. 做靈敏度分析，即各個部分的權重占比發生變化時，會給最終結果帶來什么樣的變化。

     如下的靈敏度分析代碼：

     W1 = 0.1:0.001:0.5; % 表示比重以步長為 0.001 從 0.1 變到 0.5
     W2 = 1 - W1;
     n = length(W1); % 有多少種靈敏度測試樣例
     F1 = zeros(n, 1); % 目標函數1隨着比重的變化可能有的結果
     F2 = zeros(n, 1);
     X1 = zeros(n, 1); % 決策變量可能取得值
     X2 = zeros(n, 1);
     FVAL = zeros(n, 1); % 加權組合后得值可能取得值
     %%
     %% 約束條件定義
     %%
     for i = 1:n
     	w1 = W1(i); w2 = W2(i);
     	%% 
     	%% 其他部分，如系數定義
     	%%
     	[x, fval] = ……
     	F1(i) = ……
     	……
     	FVAL(i) = fval; % 記錄下各個變量，畫圖用
     	
     end
     
     figure(1)
     plot(W1, F1, W1, F2) % 兩條線：分別是，兩個目標函數隨着 W1 比重的變化而變化
     xlabel('f_{1}的權重')
     ylabel('f_{1}和f_{2}的取值')
     legend('f_{1}','f_{2}') % 右上角對線標號
     
     figure(2)
     plot(W1, X1, W1, X2) 
     xlabel('f_{1}的權重')
     ylabel('X_{1}和X_{2}的取值')
     legend('X_{1}','X_{2}') 
     
     figure(3)
     plot(W1, FVAL)
     xlabel('f_{1}的權重')
     ylabel('綜合指標值')
     
     end