lomena_konzola_MC

Metoda Monte Carlo (MC)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import colors
import seaborn as sns
import pandas as pd
from scipy import stats
import tabulate
In [2]:
%matplotlib notebook
In [3]:
def background_gradient(self, cvals=None, cmin=None, cmax=None, cmap='viridis', **css):
    """For use with `DataFrame.style.apply` this function will apply a heatmap
    color gradient *elementwise* to the calling DataFrame

    Parameters
    ----------
    self : pd.DataFrame
        The calling DataFrame. This argument is automatically passed in by the
        `DataFrame.style.apply` method

    cvals : pd.DataFrame
        If specified this DataFrame is used to determine the color gradient

    cmin : float
        If specified, any values below this will be clipped to the bottom of
        the cmap

    cmax : float
        If specified, any values above this will be clipped to the top of
        the cmap

    cmap : colormap or str
        The colormap to use

    css : dict
        Any extra inline css key/value pars to pass to the styler

    Returns
    -------
    pd.DataFrame
        The css styles to apply

    """
    from matplotlib import cm
    from matplotlib.colors import rgb2hex
    if cvals is None:
        cvals = self.values.ravel().copy()
    else:
        assert cvals.shape == self.shape
        cvals = cvals.values.ravel().copy()
    cvals -= cmin or cvals.min()
    cvals /= cmax or cvals.max()
    cvals = cvals.clip(0, 1)
    styles = []
    for rgb in cm.get_cmap(cmap)(cvals):
        style = [
            "{}: {}".format(key.replace('_', '-'), value)
            for key, value in css.items()
        ]
        style.append("background-color: {}".format(rgb2hex(rgb)))
        styles.append('; '.join(style))
    styles = np.asarray(styles).reshape(self.shape)
    return pd.DataFrame(styles, index=self.index, columns=self.columns)

Vstupní náhodné veličiny

  • všechny veličiny mají normální rozdělení
In [4]:
rv_a1 = stats.norm(0.1, 0.01) # m
rv_a2 = stats.norm(0.2, 0.01) # m
rv_f1 = stats.norm(15000, 1500) # N
rv_f2 = stats.norm(50000, 5000) # N
rv_l1 = stats.norm(1, 0.01) # m
rv_l2 = stats.norm(2, 0.01) # m
rv_ft1 = stats.norm(103e6, 15e6) # Pa
rv_ft2 = stats.norm(103e6, 15e6) # Pa
In [5]:
# seznam náhodných veličin - vstupních a počítaných
rv_list = ['a1', 'a2', 'f1', 'f2', 'l1', 'l2', 'ft1', 'ft2', 'e1', 'e2', 'z1', 'z2', 'z']
In [6]:
ncols =  2
nrows = 4
fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(8, 6), tight_layout=True)
for i in range(nrows):
    for j in range(ncols):
        rv_name = rv_list[i * ncols + j]
        rv = locals().get('rv_' + rv_name)
        if not rv: continue
        x = np.linspace(rv.mean() - rv.std()*3, rv.mean() + rv.std()*3, 1000)
        ax[i, j].plot(x, rv.pdf(x))
        ax[i, j].set_title(rv_name)
No description has been provided for this image

Výpočet metodou MC

Pravděpodobnosti poruchy

  • $p_{f1}$ - v řezu I
  • $p_{f2}$ - v řezu II
  • $p_{f}$ - celé konstrukce (minimum($Z_1$ a $Z_2$))
In [7]:
def MC(nsim):
    # generování realizací pro jednotlivé náhodné veličiny
    a1 = rv_a1.rvs(nsim)
    a2 = rv_a2.rvs(nsim)
    f1 = rv_f1.rvs(nsim)
    f2 = rv_f2.rvs(nsim)
    l1 = rv_l1.rvs(nsim)
    l2 = rv_l2.rvs(nsim)
    ft1 = rv_ft1.rvs(nsim)
    ft2 = rv_ft2.rvs(nsim)

    # výpočet účinků od zatížení
    e1 = 6 * f1 * l1 / a1**3
    e2 = 6 * (f1 * l1 + f2 * l2) / a2**3 - f1 / a2**2

    # odolnost konstrukce
    r1 = ft1
    r2 = ft2

    # rezerva spolehlivosti
    z1 = r1 - e1
    z2 = r2 - e2
    z = np.minimum(z1, z2)
    
    # výpočet pravděpodobností poruchy
    pf_1 = np.sum(z1<0) / nsim
    pf_2 = np.sum(z2<0) / nsim
    pf = np.sum(z < 0) / nsim # porucha celé konstrukce (min)
    return a1, a2, f1, f2, l1, l2, ft1, ft2, e1, e2, r1, r2, z1, z2, z, pf_1, pf_2, pf
In [8]:
nsim = 1000
a1, a2, f1, f2, l1, l2, ft1, ft2, e1, e2, r1, r2, z1, z2, z, pf_1, pf_2, pf = MC(nsim)
In [9]:
ncols =  2
nrows = 4
fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(8, 6), tight_layout=True)
for i in range(nrows):
    for j in range(ncols):
        rv_name = rv_list[i * ncols + j]
        rv = locals().get('rv_' + rv_name)
        if not rv: continue
        x = np.linspace(rv.mean() - rv.std()*3, rv.mean() + rv.std()*3, 1000)
        ax[i, j].hist(locals()[rv_name], bins='sqrt', normed=True, label='MC', alpha=.5, ec='k')
        ax[i, j].plot(x, rv.pdf(x), label='exact')
        ax[i, j].set_title(rv_name)
ax[0, 0].legend()
No description has been provided for this image
Out[9]:

Pravděpodobnosti poruchy

In [10]:
print('pf_1 =', pf_1)
print('pf_2 =', pf_2)
print('pf =', pf)

pf_1 = 0.348
pf_2 = 0.201
pf = 0.475

Graf rezerv spolehlivosti

In [11]:
fig, ax = plt.subplots(tight_layout=True)
ax.hist(z1, bins=50, normed=True, alpha=.5, label='z1', ec='k')
ax.hist(z2, bins=50, normed=True, alpha=.5, label='z2', ec='k')
ax.legend();
No description has been provided for this image

Korelace mezi jednotlivými náhodnými veličinami

In [12]:
rv_arr = []
for i in rv_list:
    rv_arr.append(locals()[i])
rv_arr = np.array(rv_arr)
#print(tabulate.tabulate(np.corrcoef(rv_arr), floatfmt='.3f', tablefmt='grid', headers=rv_list, showindex=rv_list))
corr_tab = pd.DataFrame(np.corrcoef(rv_arr), index=rv_list, columns=rv_list)
display(corr_tab.style.apply(background_gradient, axis=None, cmap='bwr', cmin=-1, cmax=2).format('{:.3f}'))
a1 a2 f1 f2 l1 l2 ft1 ft2 e1 e2 z1 z2 z
a1 1.000 0.014 -0.036 -0.007 0.002 0.018 -0.010 -0.002 -0.910 -0.022 0.824 0.014 0.684
a2 0.014 1.000 -0.023 0.061 -0.041 -0.030 0.013 0.043 -0.032 -0.855 0.034 0.629 0.240
f1 -0.036 -0.023 1.000 -0.006 0.016 -0.019 -0.026 0.054 0.332 0.089 -0.312 -0.025 -0.270
f2 -0.007 0.061 -0.006 1.000 -0.030 0.044 0.013 0.014 -0.005 0.442 0.010 -0.300 -0.086
l1 0.002 -0.041 0.016 -0.030 1.000 -0.032 -0.020 -0.033 0.022 0.016 -0.028 -0.034 -0.029
l2 0.018 -0.030 -0.019 0.044 -0.032 1.000 0.049 0.013 -0.030 0.067 0.046 -0.038 0.049
ft1 -0.010 0.013 -0.026 0.013 -0.020 0.049 1.000 0.069 -0.026 0.001 0.415 0.047 0.311
ft2 -0.002 0.043 0.054 0.014 -0.033 0.013 0.069 1.000 0.027 -0.029 0.003 0.714 0.223
e1 -0.910 -0.032 0.332 -0.005 0.022 -0.030 -0.026 0.027 1.000 0.054 -0.920 -0.019 -0.823
e2 -0.022 -0.855 0.089 0.442 0.016 0.067 0.001 -0.029 0.054 1.000 -0.048 -0.721 -0.288
z1 0.824 0.034 -0.312 0.010 -0.028 0.046 0.415 0.003 -0.920 -0.048 1.000 0.036 0.870
z2 0.014 0.629 -0.025 -0.300 -0.034 -0.038 0.047 0.714 -0.019 -0.721 0.036 1.000 0.357
z 0.684 0.240 -0.270 -0.086 -0.029 0.049 0.311 0.223 -0.823 -0.288 0.870 0.357 1.000

Grafické znázornění vztahů mezi jednotlivými náhodnými veličinami

In [13]:
# počet možností
nvar = len(rv_list)
nvar * (nvar - 1) / 2
Out[13]:
78.0
In [14]:
ncols =  3
nrows = 26
fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(8, nrows*2), tight_layout=True)
ix, iy = 0, 0
for i in range(nvar):
    for j in range(i+1, nvar):
        rv_name_x = rv_list[i]
        rv_name_y = rv_list[j]
        x = locals()[rv_name_x]
        y = locals()[rv_name_y]
        ax[ix, iy].plot(x, y, 'b.', ms=2)
        ax[ix, iy].set_xlabel(rv_name_x)
        ax[ix, iy].set_ylabel(rv_name_y)
        ax[ix, iy].set_title('corr = {:.3f}'.format(corr_tab.loc[rv_name_x, rv_name_y]))
        ax[ix, iy].set_xticklabels([])
        ax[ix, iy].set_yticklabels([])
        iy += 1
        if iy >= ncols:
            iy = 0
            ix += 1
No description has been provided for this image

Vztah mezi rezervami spolehlivosti a ostatními náhodnými veličinami

  • červeně zvýrazněny body, pro které platí $Z<0$
Rezerva zpolehlivosti $Z_1$
In [15]:
ncols =  3
nrows = 5
fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(8, nrows*2), tight_layout=True)
mask = z1 < 0
ix, iy = 0, 0
for i in range(nvar):
    rv_name_x = 'z1'
    rv_name_y = rv_list[i]
    x = locals()[rv_name_x]
    y = locals()[rv_name_y]
    ax[ix, iy].plot(x[~mask], y[~mask], 'g.', ms=2)
    ax[ix, iy].plot(x[mask], y[mask], 'r.', ms=2)
    ax[ix, iy].set_xlabel(rv_name_x)
    ax[ix, iy].set_ylabel(rv_name_y)
    ax[ix, iy].set_title('corr = {:.3f}'.format(corr_tab.loc[rv_name_x, rv_name_y]))
    ax[ix, iy].set_xticklabels([])
    ax[ix, iy].set_yticklabels([])
    iy += 1
    if iy >= ncols:
        iy = 0
        ix += 1
No description has been provided for this image
Rezerva zpolehlivosti $Z_2$
In [16]:
ncols =  3
nrows = 5
fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(8, nrows*2), tight_layout=True)
mask = z2 < 0
ix, iy = 0, 0
for i in range(nvar):
    rv_name_x = 'z2'
    rv_name_y = rv_list[i]
    x = locals()[rv_name_x]
    y = locals()[rv_name_y]
    ax[ix, iy].plot(x[~mask], y[~mask], 'g.', ms=2)
    ax[ix, iy].plot(x[mask], y[mask], 'r.', ms=2)
    ax[ix, iy].set_xlabel(rv_name_x)
    ax[ix, iy].set_ylabel(rv_name_y)
    ax[ix, iy].set_title('corr = {:.3f}'.format(corr_tab.loc[rv_name_x, rv_name_y]))
    ax[ix, iy].set_xticklabels([])
    ax[ix, iy].set_yticklabels([])
    iy += 1
    if iy >= ncols:
        iy = 0
        ix += 1
No description has been provided for this image
Rezerva spolehlivosti $Z$
In [17]:
ncols =  3
nrows = 5
fig, ax = plt.subplots(ncols=ncols, nrows=nrows, figsize=(8, nrows*2), tight_layout=True)
mask = z < 0
ix, iy = 0, 0
for i in range(nvar):
    rv_name_x = 'z'
    rv_name_y = rv_list[i]
    x = locals()[rv_name_x]
    y = locals()[rv_name_y]
    ax[ix, iy].plot(x[~mask], y[~mask], 'g.', ms=2)
    ax[ix, iy].plot(x[mask], y[mask], 'r.', ms=2)
    ax[ix, iy].set_xlabel(rv_name_x)
    ax[ix, iy].set_ylabel(rv_name_y)
    ax[ix, iy].set_title('corr = {:.3f}'.format(corr_tab.loc[rv_name_x, rv_name_y]))
    ax[ix, iy].set_xticklabels([])
    ax[ix, iy].set_yticklabels([])
    iy += 1
    if iy >= ncols:
        iy = 0
        ix += 1
No description has been provided for this image

Odhad chyby metody MC

$$\mathrm{CoV}_{p_\mathrm{f}} = \frac{1}{\sqrt{N_\mathrm{sim} p_f}}$$

In [18]:
nsim = 100
nrun = 1000
pfs = []
for _ in range(nrun):
    a1, a2, f1, f2, l1, l2, ft1, ft2, e1, e2, r1, r2, z1, z2, z, pf_1, pf_2, pf = MC(nsim)
    pfs.append((pf_1, pf_2, pf))
pfs = np.array(pfs)
mu_pfs = pfs.mean(axis=0)
std_pfs = pfs.std(axis=0)
cov_pfs = std_pfs / mu_pfs
print('mu_pfs =', mu_pfs)
print('std_pfs =', std_pfs)
print('cov_pfs =', cov_pfs)
print('cov_eq =', 1 / np.sqrt(nsim * mu_pfs))

mu_pfs = [0.35945 0.22328 0.50029]
std_pfs = [0.04788108 0.04157934 0.04877721]
cov_pfs = [0.1332065  0.18622062 0.09749786]
cov_eq = [0.16679413 0.21162896 0.14138036]
In [19]:
fig, ax = plt.subplots(figsize=(6,3), tight_layout=True)
ax.hist(pfs[:, 0], normed=True, bins='auto',
        ec='k', alpha=.5, label='pf_1')
ax.hist(pfs[:, 1], normed=True, bins='auto',
        ec='k', alpha=.5, label='pf_2')
ax.hist(pfs[:, 2], normed=True, bins='auto',
        ec='k', alpha=.5, label='p_f')
ax.legend();
No description has been provided for this image
In [20]:
nsim = 1000
nrun = 1000
pfs = []
for _ in range(nrun):
    a1, a2, f1, f2, l1, l2, ft1, ft2, e1, e2, r1, r2, z1, z2, z, pf_1, pf_2, pf = MC(nsim)
    pfs.append((pf_1, pf_2, pf))
pfs = np.array(pfs)
mu_pfs = pfs.mean(axis=0)
std_pfs = pfs.std(axis=0)
cov_pfs = std_pfs / mu_pfs
print('mu_pfs =', mu_pfs)
print('std_pfs =', std_pfs)
print('cov_pfs =', cov_pfs)
print('cov_eq =', 1 / np.sqrt(nsim * mu_pfs))

mu_pfs = [0.360366 0.224384 0.502187]
std_pfs = [0.01496443 0.01312938 0.01589226]
cov_pfs = [0.04152563 0.05851298 0.03164611]
cov_eq = [0.05267786 0.06675811 0.04462387]
In [21]:
fig, ax = plt.subplots(figsize=(6,3), tight_layout=True)
ax.hist(pfs[:, 0], normed=True, bins='auto',
        ec='k', alpha=.5, label='pf_1')
ax.hist(pfs[:, 1], normed=True, bins='auto',
        ec='k', alpha=.5, label='pf_2')
ax.hist(pfs[:, 2], normed=True, bins='auto',
        ec='k', alpha=.5, label='p_f')
ax.legend();
No description has been provided for this image
In [ ]:
In [ ]: