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

# --- Fonctions ---
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))
        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):

    # --- Coupures limites ---
    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])))

    # --- Calcul profits ---
    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

    # --- Solution optimale ---
    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]

    # --- Équilibres ---
    ceq, veq = [], []
    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])
    check_equilibre(v1, v2)
    check_equilibre(v1, v3)
    check_equilibre(v2, v3)

    # --- Limites axes auto ---
    y_max = [v1.max(), v2.max(), v3.max()]
    y_min = [v1.min(), v2.min(), v3.min()]
    y_min_limit = min(y_min)
    y_max_limit = max(y_max) + abs(max(y_max))*0.1

    dy = 0.035 * (max(y_max) - min(y_max))

    # --- Graphe ---
    plt.figure(figsize=(9,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'Moyenne={moy}, Variance={s2}', fontsize=10)
    plt.xlabel('Teneur de coupure (%)', fontsize=12)
    plt.ylabel('Profit / t. minéralisée ($)', fontsize=12)

    # --- Labels limites ---
    labels_limites = ['C1','C2','C3']
    v_limites = [v1.max(), v2.max(), v3.max()]
    c_limites = [c1, c2, c3]
    for c, v, label in zip(c_limites, v_limites, labels_limites):
        plt.text(c, v+dy, label, fontsize=10, fontweight='bold', color='black')
        plt.vlines(c, v, v+0.8*dy, colors='black', linewidth=1.1)

    # --- Labels équilibres ---
    labels_eq = ['$c_{12}$','$c_{13}$','$c_{23}$']
    for c, v, label in zip(ceq, veq, labels_eq):
        plt.text(c, v+dy, label, fontsize=10, fontweight='bold', color='red')
        plt.vlines(c, v, v+0.8*dy, colors='red', linewidth=1.1)

    plt.xlim(0, cmax)
    plt.ylim(y_min_limit, y_max_limit)
    plt.legend(fontsize=10)
    plt.show()

    # --- Tableau récap ---
    rows = [
        ["Limite C1", f"{c1:.2f} %", f"{v1.max():.2f} $", "Mine"],
        ["Limite C2", f"{c2:.2f} %", f"{v2.max():.2f} $", "Concentrateur"],
        ["Limite C3", f"{c3:.2f} %", f"{v3.max():.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} $", "Équilibre"])
    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.3, step=0.1,description='Coût minage (m)', style=style, layout=layout),
    'y': widgets.FloatText(value=0.9, step=0.01,description='Taux récupération (y)', style=style, layout=layout),
    'p': widgets.FloatText(value=1700, step=25,description='Prix métal (p)', style=style, layout=layout),
    'k': widgets.FloatText(value=500, step=25,description='Coût fonderie (k)', style=style, layout=layout),
    'h': widgets.FloatText(value=3, step=0.1,description='Coût traitement (h)', style=style, layout=layout),
    'f': widgets.FloatText(value=20, step=0.25,description='Frais fixes (f)', style=style, layout=layout),
    'F': widgets.FloatText(value=0, step=0.25,description='Coût opportunité (F)', style=style, layout=layout),
    'M': widgets.FloatText(value=24, step=0.25,description='Capacité minage (M)', style=style, layout=layout),
    'H': widgets.FloatText(value=14, step=0.25,description='Capacité traitement (H)', style=style, layout=layout),
    'K': widgets.FloatText(value=0.22, step=0.1,description='Capacité marché (K)', style=style, layout=layout),
    'moy': widgets.FloatText(value=1.3, step=0.01,description='Moyenne (%)', style=style, layout=layout),
    's2': widgets.FloatText(value=3, step=0.05,description='Variance (%)²', style=style, layout=layout),
}

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

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
        )

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

display(ui, out)
update_plot()