Try   HackMD

適用Monte Carlo方法的fmincon解法與資料篩選

陳宥廷 DextinWed, Oct 14, 2020

tags: 筆記

來自我修最佳設計時(2018)的作業5,當時發現fmincon決定計算機率的梯度,在作業中已經提出解法。
另外針對FOSM處理非線性問題,當初作業6中使用更好的修正方案,但還是無法解決根本問題,其實在作業4的結尾我有想到可能的解決方法,但當時沒有深入研究,最近忽然想到解決方法,想和大家分享。

  1. 範例題目
  2. Monte Carlo方法簡介
  3. fmincon對Monte Carlo問題的解法
  4. 大量計算時Monte Carlo的資料篩選法

一、蒙地卡羅法 Monte Carlo Method

利用隨機撒點來模擬真實系統的不確定性

圓周率估算

以下是估算

π的簡單範例

n = 10000;
x = rand(n,1);
y = rand(n,1);
pi_e = 4*sum([x.^2 + y.^2] <= 1)/n;
% pi = 4倍面積機率

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

二、範例

原題目

minx1,x2  R+x1+x2S.T.  20x12x201204(x1+x25)2(x1x212)20x12+8x2750

可行解區 feasible region

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

原始解

function [xsol, history, flag]=solve(x0)
    lb = [0.01 0.01]; 
    ub = [10 10]; 
    options = optimoptions(@fmincon, 'OutputFcn', @outfun, 'Display', 'off'); 
    [xsol, ~, flag] = fmincon(@(x) x(1)+x(2), x0, [], [], [], [], lb, ub, @nonlcon, options);

    function stop = outfun(x,optimValues,state) 
        stop = false; 
        if state == 'iter' 
            history = [history; x optimValues.fval]; 
        end 
    end
    
    function [g, geq]=nonlcon(x)
        g1 = 20 - x(2).*(x(1).^2);
        g2 = 120 - 4*((x(1)+x(2)-5).^2) - (x(1)-x(2)-12).^2;
        g3 = x(1).^2 + 8*x(2) - 75;
        g = [g1 g2 g3];
        geq = [];
    end
end

不確定性

給定

x1
x2
的標準差為0.3,使得違反限制條件的機率小於0.0013,新的限制條件如下:

Pr[ 20x12x2>0 ]0.00130Pr[ 1204(x1+x25)2(x1x212)2>0 ]0.00130Pr[ x12+8x275>0 ]0.00130

三、fmincon的瓶頸與解法

如果只按照新的限制條件來寫,會出現如下圖的情形,無法收斂到理想的位置

function [g, geq]=nonlcon(x)
    n = 1e4; % 採樣數量
    std = 0.3;
    tolerance = 0.0013;
    X1 = normrnd(x(1), std, [n, 1]); 
    X2 = normrnd(x(2), std, [n, 1]);
    g1 = 20 - X2.*(X1.^2);
    g2 = 120 - 4*((X1+X2-5).^2) - (X1-X2-12).^2;
    g3 = X1.^2 + 8*X2 - 75;
    g = [sum(g1>0); sum(g2>0); sum(g3>0)] - tolerance*n;
    geq = [];
end

因為fmincon無法計算機率的梯度,因為機率、統計或布林函數只能計算符合或不符合的程度,其值無法有效表示當前位置,也無法計算梯度

Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →
Image Not Showing Possible Reasons
  • The image file may be corrupted
  • The server hosting the image is unavailable
  • The image path is incorrect
  • The image format is not supported
Learn More →

解法

等價性

function [g, geq]=nonlcon(x)
    ith = fix(n*0.0013) + 1; 
    g1 = 20 - X2.*(X1.^2);
    g2 = 120 - 4*((X1+X2-5).^2) - (X1-X2-12).^2;
    g3 = X1.^2 + 8*X2 - 75;
    g = sort([g1 g2 g3],'descend');
    g = g(ith,:);
    % 取第ith大的作為代表
    geq = [];
end

四、資料篩選法(有時間再寫:p)

預先處理

條件與假設

簡單函數

複雜函數

其他fmincon