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

> [name=陳宥廷 Dextin][time=Wed, Oct 14, 2020] [color=#87CEFA]
>

###### tags: `筆記`

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

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

## 一、蒙地卡羅法 Monte Carlo Method

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

### 圓周率估算

以下是估算$\pi$的簡單範例

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

 

## 二、範例

### 原題目
$$
\begin{split}
\underset{x_1, x_2\ \in \ R^+}{\min} &x_1+x_2 \\
S.T. \ \
& 20 - x_1^2x_2 \leq 0\\
& 120 - 4(x_1 +x_2 -5)^2 - (x_1 -x_2 -12)^2 \leq 0\\
& x_1^2 + 8x_2 -75 \leq 0
\end{split}
$$

### 可行解區 feasible region

 

### 原始解

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

### 不確定性

給定$x_1$和$x_2$的標準差為0.3，使得**違反限制條件的機率**小於0.0013，新的限制條件如下：

$$
\begin{split}
& Pr[\ 20 - x_1^2x_2 > 0\ ]-0.0013 \leq 0 \\
& Pr[\ 120 - 4(x_1 +x_2 -5)^2 - (x_1 -x_2 -12)^2 > 0\ ]-0.0013 \leq 0\\
& Pr[\ x_1^2 + 8x_2 -75 > 0\ ]-0.0013 \leq 0
\end{split}
$$

### 三、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無法計算機率的梯度，因為機率、統計或布林函數只能計算符合或不符合的程度，其值無法有效表示當前**位置**，也無法計算**梯度**。
 

### 解法

等價性

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