---

0.說明

• 本文是雞年新年計劃的內容之一：每週學習一個數學模型，寫一篇總結，記錄自己學到的東西和遇到的問題。
• 這些文章並不是相關模型的全面介紹，也不是從最基礎的開始，所以不一定適合數學模型的beginner，但都是些很實際的技術，希望能幫到你。
• 本文的重點是使用Monte Carlo方法求解非線性規劃。

• • 說明
  • 向大家求助的問題
  • 用Monte Carlo 方法求解如下非線性整數規劃
  • 關於randi函式和randint函式的異同
  • 如何確定選取多少點

1.向大家求助的問題

關於下面的這些問題，我迫切的需要你的幫助，哪怕是一些網址和資訊，所以，感謝你給我發郵件[email protected]或直接加我QQ好友：117829132（請一定註明“幫助你”）或直接在評論中留言。

1.Monte Carlo方法求解非線性規劃，取多少個點是至關重要的，歡迎你告訴我更多實用的、確定點數的方法。

2.一個線性規劃模型記為M，限制其中某些（或全部）決策變數只能取整數後得到一個整數線性規劃，記為M’。即使M有最優解，M’也有可能甚至沒有可行解，為什麼？

2. 用Monte Carlo 方法求解如下非線性整數規劃

minz=x21+x22+3x23+4x24+2x25−8x1−2x2−3x3−x4−2x5 $$min \quad z = x_1^2+x_2^2+3x_3^2+4x_4^2+2x_5^2-8x_1-2x_2-3x_3-x_4-2x_5$$
$$s.t.\begin{cases} \sum_{i=1}^5x_i \le 400 \\ x_1+2x_2+2x_3+x_4+6x_5 \le 800 \\ 2x_1+x_2+6x_3 \le 200 \\ x_3+x_4+5x_5 \le 200\\ 0 \le x_i \le 99, \; i = 1, \cdots , 5 \end{cases}$$

使用MATLAB編寫程式求解，需要編寫兩個程式：

計算給定點的目標函式值和約束條件值的程式如下：

function [ f, g ] = ex2_6( x )
f = [1 1 3 4 2] * (x.^2)' - [8 2 3 1 2] * x';
g = [ sum(x) - 400
    [1 2 2 1 6] * x' - 800
    [2 1 6 0 0] * x' - 200
    [0 0 1 1 5] * x' - 200];

隨機選取點進行驗證和統計的程式如下：

fix(clock)  %顯示當前時刻，也就是程式開始執行的時刻
x00=[];p00=[];t=[];%初始化三個變數，用於存放最優解（x00）、最優值（p00）和每個單次（取10^6個點）進行計算的時間，用於後續的分析
for k=1:3000    %為了測試效能和保證準確性，進行3000個單次的迴圈
rand( 'state', sum(clock)); %初始化隨機數的種子
tic;        %開始計時
p0 = 0; x0 = 0; %初始化這個單次的最優值和最優解
for i = 1:10^6  %共選取10^6個點
    x = randi(99, 1, 5);    %隨時生成一個1行5列的x，每個元素是1至99之間的整數
    [f, g] = ex2_6(x);      %計算目標函式和約束條件的值
    if all(g <= 0)          %如果約束條件全部滿足
        if f>p0
            p0 = f; x0 = x; %如果該值處的目標函式值比當前的最優值還大，那麼更新最優值和最優解
        end
    end
end
%x0, p0
x00(k, :) = x0;     %記錄每個單次的最優解
p00(k, :) = p0;     %記錄每個單次的最優值
t(k)=toc;           %記錄每個單次的執行時間
end
x00
p00
fix(clock)          %顯示全部程式結束的時刻
t

上述程式共執行3000個單次，每個單次選擇10^6個點進行驗證。3000個單次共耗時36815秒、合計613.5793分鐘、10.2263小時，比事前預估的時間少30幾分鐘。

用Monte Carlo方法找到的最優解是[41 99 3 99 17]，最優值是，50623

該問題也可以用Lingo精確求解，準確的最優解是：[50 99 0 99 20]，最優值是51568。
二者之間的差約為1.83%，對於很多問題而言，精度已經足夠了。

2.關於randi函式和randint函式的異同

第一版的程式中，我使用的是randint函式來生成隨機數，每個單次執行時間大約1400多秒，合計24分鐘左右。每個單次選取10^6個點進行驗證。

隨後，按照MATLAB的提示，我換成randi函式，則每個單次的執行時間劇減為13至14秒之間。

開始我以為是對於程式的其他修改導致的，於是，我做了兩個版本的程式，這兩個版本之間只有該函式不同。嘗試執行後randi版本的時間仍然是14秒左右。

為了保險起見，我還在重啟電腦前後測試了randi版本的執行時間，結果是沒有變化，也就是說，可以排除之前執行的程式可能殘留的初值導致的執行時間的減少。

所以，我的結論是，使用randi函式生成隨機整數，執行時間是使用較早版本MATLAB中的randint的時間1.1%，可以大幅度的提高模擬效率！！

3. 如何確定選取多少點

用Monte Carlo方法求解非線性規劃時，一個非常重要的任務就是要確定選取多少點，點數過多則計算時間過長；點數太少則結果的準確性就會較差，如何取點是非常重要的。上述問題是純整數規劃，每個x可以取0至99共100個點，所以，全部驗證所有可能性需要(100)^5=10^10個點，需要耗費大量的時間。

本文給出的，是驗證10^6個點，下面粗略的說明，為什麼驗證10^6個點就足夠了？（不失一般性，可以建設最優值點不是孤立的奇點）

假設每個點處的函式值落在高值區的可能性為0.01和0.00001，則當驗證10^6個點後，至少有一個點落在高值區的概率為

$$1-0.99^{1000000} \approx 0.99 \dots 99（100多位）\\ 1-0.99999^{1000000} \approx 0.999954602$$
可見，我們驗證10^6個點就可以在高值區找到一個點的可能性還是非常高的。

需要注意的是，上述敘述中使用的是“高值區”，而不是最優點，這二者是不同的，請仔細區分。

未完待續