```python
import math
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from scipy.optimize import minimize
from sklearn import metrics
import time
import pandas as pd
```

###  a) The Chierici’s relative permeability model 

**Eq. 5**
\begin{equation}
\kappa_{r w}=\kappa_{w}^{0} e^{-B\left(\frac{S_{w}-S_{w 0}}{1-S_{w}-S_{o 0}}\right)^{-M}}
\quad \kappa_{r o}=\kappa_{o}^{0} e^{-A\left(\frac{S_{w}-S_{w 0}}{1-S_{w}-S_{o 0}}\right)^{L}}
\end{equation}


```python
#water relative permeability at the residual oil saturation /// krwe in code water rel perm at water curve end point
kw0=0.375
#oil relative permeability at the connate water saturation /// kroe in code oil rel perm at oil curve end point
ko0=0.95
#𝑆w0  - connate water saturation /// residual_w = residual water saturation, Swirr (irreducible)
sw0=0.2
#𝑆𝑜0 - residual oil saturation/// residual_o =, Sor
so0=0.26
    
```


```python

```



**Eq. 5**
\begin{equation}
\kappa_{r w}=\kappa_{w}^{0} e^{-B\left(\frac{S_{w}-S_{w 0}}{1-S_{w}-S_{o 0}}\right)^{-M}}
\quad \kappa_{r o}=\kappa_{o}^{0} e^{-A\left(\frac{S_{w}-S_{w 0}}{1-S_{w}-S_{o 0}}\right)^{L}}
\end{equation}


```python
from scipy import optimize
#initial data
kw0=0.375
ko0=0.95
sw0=0.2
so0=0.26

filepath = '/Users/mikhailgasanov/Desktop/Geostatistics and Reservoir Simulation/Homework/dataset1.xlsx'
qi = pd.read_excel(filepath)
```


```python
# add random table for predictions
qi['Krw_predicted'] = qi['Krw'] - 0.1
```


```python
np.sqrt(np.mean(qi['Krw']-qi['Krw_predicted'])**2)
```




    0.10000000000000003




```python
qi
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Date</th>
      <th>Sw</th>
      <th>Krw</th>
      <th>Kro</th>
      <th>Krw_predicted</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1</td>
      <td>0.20</td>
      <td>0.000</td>
      <td>0.950</td>
      <td>-0.100</td>
    </tr>
    <tr>
      <th>1</th>
      <td>2</td>
      <td>0.22</td>
      <td>0.000</td>
      <td>0.941</td>
      <td>-0.100</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3</td>
      <td>0.24</td>
      <td>0.000</td>
      <td>0.917</td>
      <td>-0.100</td>
    </tr>
    <tr>
      <th>3</th>
      <td>4</td>
      <td>0.26</td>
      <td>0.000</td>
      <td>0.879</td>
      <td>-0.100</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5</td>
      <td>0.28</td>
      <td>0.000</td>
      <td>0.828</td>
      <td>-0.100</td>
    </tr>
    <tr>
      <th>5</th>
      <td>6</td>
      <td>0.30</td>
      <td>0.000</td>
      <td>0.766</td>
      <td>-0.100</td>
    </tr>
    <tr>
      <th>6</th>
      <td>7</td>
      <td>0.32</td>
      <td>0.001</td>
      <td>0.696</td>
      <td>-0.099</td>
    </tr>
    <tr>
      <th>7</th>
      <td>8</td>
      <td>0.34</td>
      <td>0.002</td>
      <td>0.621</td>
      <td>-0.098</td>
    </tr>
    <tr>
      <th>8</th>
      <td>9</td>
      <td>0.36</td>
      <td>0.002</td>
      <td>0.545</td>
      <td>-0.098</td>
    </tr>
    <tr>
      <th>9</th>
      <td>10</td>
      <td>0.38</td>
      <td>0.004</td>
      <td>0.470</td>
      <td>-0.096</td>
    </tr>
    <tr>
      <th>10</th>
      <td>11</td>
      <td>0.40</td>
      <td>0.005</td>
      <td>0.399</td>
      <td>-0.095</td>
    </tr>
    <tr>
      <th>11</th>
      <td>12</td>
      <td>0.42</td>
      <td>0.007</td>
      <td>0.335</td>
      <td>-0.093</td>
    </tr>
    <tr>
      <th>12</th>
      <td>13</td>
      <td>0.44</td>
      <td>0.010</td>
      <td>0.277</td>
      <td>-0.090</td>
    </tr>
    <tr>
      <th>13</th>
      <td>14</td>
      <td>0.46</td>
      <td>0.014</td>
      <td>0.227</td>
      <td>-0.086</td>
    </tr>
    <tr>
      <th>14</th>
      <td>15</td>
      <td>0.48</td>
      <td>0.018</td>
      <td>0.184</td>
      <td>-0.082</td>
    </tr>
    <tr>
      <th>15</th>
      <td>16</td>
      <td>0.50</td>
      <td>0.024</td>
      <td>0.147</td>
      <td>-0.076</td>
    </tr>
    <tr>
      <th>16</th>
      <td>17</td>
      <td>0.52</td>
      <td>0.031</td>
      <td>0.116</td>
      <td>-0.069</td>
    </tr>
    <tr>
      <th>17</th>
      <td>18</td>
      <td>0.54</td>
      <td>0.040</td>
      <td>0.091</td>
      <td>-0.060</td>
    </tr>
    <tr>
      <th>18</th>
      <td>19</td>
      <td>0.56</td>
      <td>0.051</td>
      <td>0.070</td>
      <td>-0.049</td>
    </tr>
    <tr>
      <th>19</th>
      <td>20</td>
      <td>0.58</td>
      <td>0.064</td>
      <td>0.053</td>
      <td>-0.036</td>
    </tr>
    <tr>
      <th>20</th>
      <td>21</td>
      <td>0.60</td>
      <td>0.080</td>
      <td>0.040</td>
      <td>-0.020</td>
    </tr>
    <tr>
      <th>21</th>
      <td>22</td>
      <td>0.62</td>
      <td>0.100</td>
      <td>0.029</td>
      <td>0.000</td>
    </tr>
    <tr>
      <th>22</th>
      <td>23</td>
      <td>0.64</td>
      <td>0.123</td>
      <td>0.021</td>
      <td>0.023</td>
    </tr>
    <tr>
      <th>23</th>
      <td>24</td>
      <td>0.66</td>
      <td>0.151</td>
      <td>0.014</td>
      <td>0.051</td>
    </tr>
    <tr>
      <th>24</th>
      <td>25</td>
      <td>0.68</td>
      <td>0.184</td>
      <td>0.009</td>
      <td>0.084</td>
    </tr>
    <tr>
      <th>25</th>
      <td>26</td>
      <td>0.70</td>
      <td>0.222</td>
      <td>0.006</td>
      <td>0.122</td>
    </tr>
    <tr>
      <th>26</th>
      <td>27</td>
      <td>0.72</td>
      <td>0.267</td>
      <td>0.003</td>
      <td>0.167</td>
    </tr>
    <tr>
      <th>27</th>
      <td>28</td>
      <td>0.74</td>
      <td>0.317</td>
      <td>0.002</td>
      <td>0.217</td>
    </tr>
    <tr>
      <th>28</th>
      <td>29</td>
      <td>0.76</td>
      <td>0.375</td>
      <td>0.001</td>
      <td>0.275</td>
    </tr>
  </tbody>
</table>
</div>




```python
class Project():
    def __init__(self, df):
        self.df = df
        self.kw0=0.375
        self.ko0=0.95
        self.sw0=0.2
        self.so0=0.26
    
    def k_rw(self, parameters):
        #compute k_rw for each sw in dataframe
        k_rw_predicted = []
        B, M = parameters

        k_rw = self.kw0
        for i in range(len(test.df)):
            if i==0:
                _part = (qi.loc[i, 'Sw']-self.sw0)/(1-qi.loc[i, 'Sw']-self.so0)
                k_rw_new = self.kw0 * np.exp((-B*_part)**(-M))
            else: 
                _part = (qi.loc[i, 'Sw']-qi.loc[i-1, 'Sw'])/(1-qi.loc[i, 'Sw']-self.so0)
                k_rw_new = k_rw * np.exp((-B*_part)**(-M))
            k_rw_predicted.append(k_rw_new)
            k_rw = k_rw_new

        rmse = np.sqrt(np.mean(self.df['Krw']-k_rw_predicted)**2)
        print(rmse)
        return rmse
    def optimization(self, B, M):
        x0 = np.array([B,M])
        bnds = ((1.86,2.96), (0.64, 0.97))
        res = optimize.minimize(self.k_rw, x0, bounds=bnds, \
        method='COBYLA')
        return res
```


```python
# B=2.83
# M=0.64
test = Project(qi)
B, M = 2.83, 0.64
res = test.optimization(B,M)
res
```

    inf
    inf
    inf


    /Users/mikhailgasanov/miniconda3/lib/python3.8/site-packages/scipy/optimize/_minimize.py:541: RuntimeWarning: Method COBYLA cannot handle bounds.
      warn('Method %s cannot handle bounds.' % method,
    <ipython-input-20-e2359d801e01>:18: RuntimeWarning: divide by zero encountered in double_scalars
      k_rw_new = self.kw0 * np.exp((-B*_part)**(-M))
    <ipython-input-20-e2359d801e01>:21: RuntimeWarning: invalid value encountered in double_scalars
      k_rw_new = k_rw * np.exp((-B*_part)**(-M))
    <ipython-input-20-e2359d801e01>:20: RuntimeWarning: divide by zero encountered in double_scalars
      _part = (qi.loc[i, 'Sw']-qi.loc[i-1, 'Sw'])/(1-qi.loc[i, 'Sw']-self.so0)





         fun: inf
       maxcv: 0.0
     message: 'NaN result encountered.'
        nfev: 3
      status: 5
     success: False
           x: array([2.83, 0.64])




```python
#Trying to compute without class
k_rw_predicted = []
B, M = 2.83, 0.64
k_rw = test.kw0
for i in range(len(test.df)):
    sw = qi.loc[i, 'Sw']
    if i==0:
        k_rw_new = test.kw0 * np.exp((-B*((qi.loc[i, 'Sw']-test.sw0)/(1-qi.loc[i, 'Sw']-test.so0)))**(-M))
    else: 
        _part = (qi.loc[i, 'Sw']-qi.loc[i-1, 'Sw'])/(1-qi.loc[i, 'Sw']-test.so0)
#         print(np.exp(-B*_part**(-M)))
        k_rw_new = k_rw * np.exp(-B*_part**(-M))
    k_rw_predicted.append(k_rw_new)
    k_rw = k_rw_new
#     print(k_rw)
rmse = np.sqrt(np.mean(test.df['Krw']-k_rw_predicted)**2)
```

    <ipython-input-23-32a231bc5eaf>:7: RuntimeWarning: divide by zero encountered in double_scalars
      k_rw_new = test.kw0 * np.exp((-B*((qi.loc[i, 'Sw']-test.sw0)/(1-qi.loc[i, 'Sw']-test.so0)))**(-M))
    <ipython-input-23-32a231bc5eaf>:9: RuntimeWarning: divide by zero encountered in double_scalars
      _part = (qi.loc[i, 'Sw']-qi.loc[i-1, 'Sw'])/(1-qi.loc[i, 'Sw']-test.so0)
    <ipython-input-23-32a231bc5eaf>:11: RuntimeWarning: invalid value encountered in double_scalars
      k_rw_new = k_rw * np.exp(-B*_part**(-M))



```python
rmse
```




    inf




```python
k_rw * np.exp((-B*((qi.loc[i, 'Sw']-qi.loc[i-1, 'Sw'])/(1-qi.loc[i, 'Sw']-test.so0)))**(-M))
```




    nan



5c084edee1e9c4e69c964be2843946558c458cc9
ipynb
Jupyter Notebook
Project_GeoStatistics_RP models.ipynb
ABAZH666/project_uq_and_sa
0b47bbcb52748935cd8baecdd2129148ddfd8f19
CC0-1.0
Project_GeoStatistics_RP models.ipynb
ABAZH666/project_uq_and_sa
0b47bbcb52748935cd8baecdd2129148ddfd8f19
CC0-1.0
Project_GeoStatistics_RP models.ipynb
ABAZH666/project_uq_and_sa
0b47bbcb52748935cd8baecdd2129148ddfd8f19
CC0-1.0
Qwen/Qwen-72B
1. YES
2. YES