Atelier 1 - Feuillet de calcul — Teneur de coupure optimale

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

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

def lane(m, y, p, k, h, f, F, M, H, K, moy, s2, x_expand=1, y_expand=5):
    cmax = 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 = c_opt - x_expand, c_opt + x_expand
    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))

# Widgets
style = {'description_width': '140px'}
layout = widgets.Layout(width='210px')

params = {
    'm': widgets.FloatText(value=1.32, description='Coût minage (m)', style=style, layout=layout),
    'y': widgets.FloatText(value=0.87, description='Taux récupération (y)', style=style, layout=layout),
    'p': widgets.FloatText(value=600, description='Prix métal (p)', style=style, layout=layout),
    'k': widgets.FloatText(value=0, description='Coût fonderie (k)', style=style, layout=layout),
    'h': widgets.FloatText(value=3.41, description='Coût traitement (h)', style=style, layout=layout),
    'f': widgets.FloatText(value=11.9, description='Frais fixes (f)', style=style, layout=layout),
    'F': widgets.FloatText(value=0, description='Coût opportunité (F)', style=style, layout=layout),
    'M': widgets.FloatText(value=12, description='Capacité minage (M)', style=style, layout=layout),
    'H': widgets.FloatText(value=3.9, description='Capacité traitement (H)', style=style, layout=layout),
    'K': widgets.FloatText(value=0.085, description='Capacité marché (K)', style=style, layout=layout),
    'moy': widgets.FloatText(value=1, description='Moyenne (%)', style=style, layout=layout),
    's2': widgets.FloatText(value=4, description='Variance (%)²', style=style, layout=layout),
    'x_expand': widgets.FloatText(value=1, description='Élargissement X', style=style, layout=layout),
    'y_expand': widgets.FloatText(value=1, description='Élargissement Y', style=style, layout=layout)
}

ui = widgets.VBox([
    widgets.HBox([params['m'], params['y'], params['p']]),
    widgets.HBox([params['k'], params['h'], params['f']]),
    widgets.HBox([params['F'], params['M'], params['H']]),
    widgets.HBox([params['K'], params['moy'], params['s2']]),
    widgets.HBox([params['x_expand'], params['y_expand']])
])


out = widgets.Output()

def update_plot(change=None):
    with out:
        out.clear_output(wait=True)
        lane(
            m=params['m'].value, y=params['y'].value, p=params['p'].value, k=params['k'].value,
            h=params['h'].value, f=params['f'].value, F=params['F'].value, M=params['M'].value,
            H=params['H'].value, K=params['K'].value, moy=params['moy'].value, s2=params['s2'].value,
            x_expand=params['x_expand'].value, y_expand=params['y_expand'].value
        )

for widget in params.values():
    widget.observe(update_plot, names='value')

display(ui, out)
update_plot()