Atelier 3 - - Analyse de sensibilité

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import pandas as pd
import ipywidgets as widgets
from IPython.display import display, clear_output

# Fonction de calcul des réserves
def reserve(moy, s2, c, itype=1):
    c = np.array(c)
    m = moy
    if itype == 1:
        b = np.sqrt(np.log(s2 / m**2 + 1))
        u = np.log(m) - 0.5 * b**2
        n1 = (np.log(m / c) / b) - b / 2
        n2 = n1 + b
        tc = norm.cdf(n1)
        qc = m * norm.cdf(n2)
        mc = qc / tc
    else:
        s = np.sqrt(s2)
        tc = norm.cdf((m - c) / s)
        qc = m * tc + s * norm.pdf((c - m) / s)
        mc = qc / tc
    return tc, qc, mc

# Fonction principale
def lane(**kwargs):
    m, y, p, k, h, f, F, M, H, K, moy, s2 = [kwargs[k] for k in ['m', 'y', 'p', 'k', 'h', 'f', 'F', 'M', 'H', 'K', 'moy', 's2']]
    x_expand, y_expand = kwargs.get('x_expand', 1), kwargs.get('y_expand', 5)

    cmax = 4 #moy + 0.5 * np.sqrt(s2)
    cc = np.linspace(0.00001, cmax, 5000)

    c1 = h / (y * (p - k)) * 100
    c2 = (h + (f + F) / H) / ((p - k) * y) * 100
    c3 = h / ((p - k) - (f + F) / K) / y * 100

    cc = np.unique(np.concatenate((cc, [c1, c2, c3])))

    tc, qc, mc = reserve(moy, s2, cc, itype=1)
    tc_mc_y_pct = tc * mc * y / 100

    v1 = (p - k) * tc_mc_y_pct - tc * h - m - (f + F) / M
    v2 = (p - k) * tc_mc_y_pct - tc * h - m - (f + F) * tc / H
    v3 = (p - k) * tc_mc_y_pct - tc * h - m - (f + F) * tc_mc_y_pct / K

    v_min = np.minimum(np.minimum(v1, v2), v3)
    idx_opt = np.argmax(v_min)
    c_opt = cc[idx_opt]
    profit_opt = v_min[idx_opt]

    ceq, veq, ind = [], [], []

    def check_equilibre(vA, vB):
        diff = np.abs(vA - vB)
        i_val = np.argmin(diff)
        seuil = 10 * np.max(np.abs(np.diff(v1)))
        if diff[i_val] < seuil:
            ceq.append(cc[i_val])
            veq.append(vA[i_val])
            ind.append(i_val)

    check_equilibre(v1, v2)
    check_equilibre(v1, v3)
    check_equilibre(v2, v3)

    limites = [c1, c2, c3]
    t = np.array(limites + ceq)
    i = np.argmin(np.abs(c_opt - t))
    c_opt = t[i]

    labels_limites = ['Mine Limite', 'Concentrateur Limite', 'Marché Limite']
    labels_equilibre_possibles = ['Équilibre Mine-Concentrateur', 'Équilibre Mine-Marché', 'Équilibre Concentrateur-Marché']
    labels = labels_limites + labels_equilibre_possibles[:len(ceq)]
    nature = labels[i]

    xmin, xmax = 0, 3
    ymin, ymax = profit_opt - y_expand, profit_opt + y_expand
    dy = 0.035 * (np.max([v1.max(), v2.max(), v3.max()]) - np.min([v1.min(), v2.min(), v3.min()]))

    plt.figure(figsize=(8, 5))
    plt.plot(cc, v1, '-', label='Mine')
    plt.plot(cc, v2, '--', label='Concentrateur')
    plt.plot(cc, v3, ':', label='Marché')
    plt.grid(True)

    min_curve = np.minimum(np.minimum(v1, v2), v3)
    plt.fill_between(cc, -100, min_curve, color='grey', alpha=0.3, hatch='///')

    plt.title(f'y={y:.2f}; p={p}; k={k}; h={h}; m={m}; f={f}; F={F}; M={M}; H={H}; K={K}; Moyenne={moy}, Variance={s2}', fontsize=10)
    plt.xlabel('Teneur de coupure (%)', fontsize=12)
    plt.ylabel('Profit / t. minéralisée ($)', fontsize=12)

    for v, label in zip([v1, v2, v3], ['$c_1$', '$c_2$', '$c_3$']):
        idx = np.argmax(v)
        if xmin <= cc[idx] <= xmax and ymin <= v[idx] <= ymax:
            plt.text(cc[idx], v[idx] + dy, label, fontsize=10, fontweight='bold', color='black')
            plt.vlines(cc[idx], v[idx], v[idx] + 0.8 * dy, colors='black', linewidth=1.1)

    pairs = [(v1, v2, '$c_{12}$'), (v1, v3, '$c_{13}$'), (v2, v3, '$c_{23}$')]
    for vA, vB, label in pairs:
        diff = np.abs(vA - vB)
        idx = np.argmin(diff)
        seuil = 0.1 * (np.max([vA.max(), vB.max()]) - np.min([vA.min(), vB.min()]))
        if diff[idx] < seuil and xmin <= cc[idx] <= xmax and ymin <= vA[idx] <= ymax:
            plt.text(cc[idx], vA[idx] + dy, label, fontsize=10, fontweight='bold', color='red')
            plt.vlines(cc[idx], vA[idx], vA[idx] + 0.8 * dy, colors='red', linewidth=1.1)

    plt.xlim(xmin, xmax)
    plt.ylim(ymin, ymax)
    plt.legend(fontsize=10)
    plt.show()

    rows = [
        ["Optimale", f"{c_opt:.2f} %", f"{profit_opt:.2f} $", nature],
        ["Limite C1", f"{c1:.2f} %", f"{v1[np.argmax(v1)]:.2f} $", "Mine"],
        ["Limite C2", f"{c2:.2f} %", f"{v2[np.argmax(v2)]:.2f} $", "Concentrateur"],
        ["Limite C3", f"{c3:.2f} %", f"{v3[np.argmax(v3)]:.2f} $", "Marché"]
    ]

    for k, (i, j) in enumerate([(0, 1), (0, 2), (1, 2)]):
        if k < len(ceq):
            rows.append([
                f"Équilibre C{i+1}-C{j+1}",
                f"{ceq[k]:.2f} %",
                f"{veq[k]:.2f} $",
                labels_equilibre_possibles[k]
            ])

    df = pd.DataFrame(rows, columns=["Type de coupure", "Teneur", "Profit associé", "Remarque"])
    print(df.to_string(index=False))

# Valeurs par défaut
default_params = {
    'm': 1.32, 'y': 0.87, 'p': 600, 'k': 0, 'h': 3.41, 'f': 11.9, 'F': 0,
    'M': 12, 'H': 3.9, 'K': 0.085, 'moy': 1, 's2': 4, 'x_expand': 1, 'y_expand': 1
}

# Widgets
param_names = ['m', 'y', 'p', 'k', 'h', 'f', 'F', 'M', 'H', 'K', 'moy', 's2']
param_labels = {
    'm': 'Coût minage (m)', 'y': 'Taux récupération (y)', 'p': 'Prix métal (p)',
    'k': 'Coût fonderie (k)', 'h': 'Coût traitement (h)', 'f': 'Frais fixes (f)',
    'F': 'Coût opportunité (F)', 'M': 'Capacité minage (M)', 'H': 'Capacité traitement (H)',
    'K': 'Capacité marché (K)', 'moy': 'Moyenne (%)', 's2': 'Variance (%)²'
}

selector = widgets.Dropdown(options=param_names, description='Paramètre à tester:')
value_slider = widgets.FloatSlider(min=0.1, max=1000, step=0.1, description='Valeur:')

out = widgets.Output()

def update_slider_range(change):
    param = selector.value
    val = default_params[param]
    value_slider.description = param_labels[param]
    value_slider.value = val
    if param in ['moy']:
        value_slider.min = 0.01
        value_slider.max = 2
        value_slider.step = 0.01
    elif param in ['s2']:
        value_slider.min = 0.01
        value_slider.max = 10
        value_slider.step = 0.01
    elif param in ['y']:
        value_slider.min = 0.2
        value_slider.max = 1
        value_slider.step = 0.01
    else:
        value_slider.min = 0
        value_slider.max = 1000
        value_slider.step = 1

def update_plot(change=None):
    param = selector.value
    val = value_slider.value
    params = default_params.copy()
    params[param] = val
    with out:
        clear_output(wait=True)
        lane(**params)

selector.observe(update_slider_range, names='value')
value_slider.observe(update_plot, names='value')

update_slider_range(None)
update_plot()

display(widgets.VBox([selector, value_slider]), out)