Ottimizzazione Multi-Obiettivo: Come Bilanciare Profitto, Rischio e Sostenibilità con Python

1. Comprensione del problema

Abbiamo tre progetti di investimento:

• A (Solare): rendimento 10%, rischio 40, [math]CO_2[/math] risparmiata 500 kg/anno (per migliaia di euro investiti)
• B (Eolico): rendimento 6%, rischio 20, [math]CO_2[/math] risparmiata 800 kg/anno
• C (Idrogeno): rendimento 4%, rischio 60, [math]CO_2[/math] risparmiata 1500 kg/anno

Capitale totale disponibile: 100 000 €.

Obiettivi (in ordine di importanza deciso dallo stakeholder):

• Massimizzare il profitto [math]f_1[/math]
• Minimizzare il rischio [math]f_2[/math] (ma qui lo trasformiamo in un vincolo)
• Massimizzare il risparmio di [math]CO_2[/math] [math]f_3[/math] (anch’esso diventa un vincolo)

Lo stakeholder fissa due soglie:

1. Il rischio medio del portafoglio non deve superare 35 (cioè la media ponderata dei rischi [math]\le 35[/math])
2. Il risparmio di [math]CO_2[/math] deve essere almeno 80 000 kg

2. Modellazione matematica

Indichiamo con [math]x_A, x_B, x_C[/math] gli importi investiti (in migliaia di euro) nei rispettivi progetti.

Il vincolo di bilancio è:
[math]x_A + x_B + x_C = 100 \text{ (in migliaia di €)}[/math]

Le funzioni obiettivo (da massimizzare/minimizzare) sono:

[math]\displaystyle \begin{aligned}
f_1 &= 0,10x_A + 0,06x_B + 0,04x_C \\
&{} \quad \text{(profitto in migliaia di €)} \\
f_2 &= 40x_A + 20x_B + 60x_C \\
&{} \quad \text{(punteggio di rischio totale)} \\
f_3 &= 500x_A + 800x_B + 1500x_C \\
&{} \quad \text{(CO}_2 \text{ risparmiata in kg)}
\end{aligned}[/math]

3. Applicazione del metodo [math]\epsilon[/math]-constraint

Scegliamo di massimizzare il profitto [math]f_1[/math] e trasformiamo gli altri due obiettivi in vincoli:

Rischio medio [math]\le 35 \implies \frac{40x_A + 20x_B + 60x_C}{x_A + x_B + x_C} \le 35[/math]

Dato che il denominatore è 100, otteniamo:
[math]40x_A + 20x_B + 60x_C \le 35 \times 100 = 3500[/math]

[math]CO_2[/math] risparmiata [math]\ge 80\,000[/math] kg:
[math]500x_A + 800x_B + 1500x_C \ge 80\,000[/math]

Aggiungiamo i vincoli di non negatività: [math]x_A, x_B, x_C \ge 0[/math].

4. Riduzione a due variabili

Usiamo l’uguaglianza [math]x_C = 100 – x_A – x_B[/math] (con [math]x_A + x_B \le 100[/math]) e sostituiamo in tutti i vincoli.

Vincolo di rischio:

[math]\displaystyle \begin{aligned}
&40x_A + 20x_B + 60(100 – x_A – x_B) \le 3500 \\
&40x_A + 20x_B + 6000 – 60x_A – 60x_B \le 3500 \\
&-20x_A – 40x_B \le -2500 \implies 20x_A + 40x_B \ge 2500 \\
&\text{Dividendo per 20:} \\
&x_A + 2x_B \ge 125 \quad (1)
\end{aligned}[/math]

Vincolo di [math]CO_2[/math]:

[math]\displaystyle \begin{aligned}
&500x_A + 800x_B + 1500(100 – x_A – x_B) \ge 80\,000 \\
&500x_A + 800x_B + 150\,000 – 1500x_A – 1500x_B \ge 80\,000 \\
&-1000x_A – 700x_B \ge -70\,000 \implies 1000x_A + 700x_B \le 70\,000 \\
&\text{Dividendo per 100:} \\
&10x_A + 7x_B \le 700 \quad (2)
\end{aligned}[/math]

Vincolo di non negatività di [math]x_C[/math]:
[math]x_A + x_B \le 100 \quad (3)[/math]
Inoltre [math]x_A \ge 0, x_B \ge 0[/math].

5. Funzione obiettivo da massimizzare

Esprimiamo [math]f_1[/math] in funzione di [math]x_A, x_B[/math]:

[math]\displaystyle \begin{aligned}
f_1 &= 0,10x_A + 0,06x_B + 0,04(100 – x_A – x_B) \\
&= 4 + 0,06x_A + 0,02x_B
\end{aligned}[/math]

Massimizzare [math]f_1[/math] equivale a massimizzare [math]0,06x_A + 0,02x_B[/math], cioè [math]3x_A + x_B[/math] (a meno di un fattore costante).

6. Determinazione della regione ammissibile

Disegniamo (idealmente) le rette nel piano [math](x_A, x_B)[/math]:

• (1) [math]x_A + 2x_B = 125[/math] (passa per (0; 62,5) e (125; 0))
• (2) [math]10x_A + 7x_B = 700[/math] (passa per (0; 100) e (70; 0))
• (3) [math]x_A + x_B = 100[/math] (passa per (0; 100) e (100; 0))

La regione ammissibile è l’insieme dei punti che soddisfano:

[math]\displaystyle \begin{cases} x_A + 2x_B \ge 125 \\ 10x_A + 7x_B \le 700 \\ x_A + x_B \le 100 \\ x_A \ge 0, x_B \ge 0 \end{cases}[/math]

Intersezione tra (1) e (2):

[math]\displaystyle \begin{aligned}
&\begin{cases} x_A + 2x_B = 125 \\ 10x_A + 7x_B = 700 \end{cases} \\
&x_A = 125 – 2x_B \implies 10(125 – 2x_B) + 7x_B = 700 \\
&1250 – 20x_B + 7x_B = 700 \implies 13x_B = 550 \\
&x_B = 550/13 \approx 42,3077; \quad x_A \approx 40,3846 \\
&\text{Punto } P_1 (40,38; 42,31) \text{. Verifica (3): } 82,69 \le 100 \text{ (OK).}
\end{aligned}[/math]

Intersezione tra (2) e (3):

[math]\displaystyle \begin{aligned}
&\begin{cases} 10x_A + 7x_B = 700 \\ x_A + x_B = 100 \end{cases} \\
&10(100 – x_B) + 7x_B = 700 \implies 3x_B = 300 \\
&x_B = 100, \quad x_A = 0 \\
&\text{Punto } P_2 (0; 100) \text{. Verifica (1): } 200 \ge 125 \text{ (OK).}
\end{aligned}[/math]

Intersezione tra (1) e l’asse [math]x_A = 0[/math]:

[math]\displaystyle \begin{aligned}
&x_B = 62,5; \quad x_A = 0 \\
&\text{Punto } P_3 (0; 62,5) \text{. Verifica (2) e (3): OK.}
\end{aligned}[/math]

7. Massimizzazione della funzione obiettivo

Valutiamo [math]0,06x_A + 0,02x_B[/math] nei vertici:

• [math]P_3(0; 62,5): 0,06(0) + 0,02(62,5) = 1,25[/math]
• [math]P_1(40,38; 42,31): 0,06(40,38) + 0,02(42,31) = 3,269[/math]
• [math]P_2(0; 100): 0,06(0) + 0,02(100) = 2,00[/math]

Il massimo si raggiunge in [math]P_1[/math].

8. Soluzione ottima

Investimenti (in migliaia di euro):

[math]\displaystyle \begin{aligned}
x_A &= 40,38 \quad (\approx 40\,380 \text{ €}) \\
x_B &= 42,31 \quad (\approx 42\,310 \text{ €}) \\
x_C &= 100 – 40,38 – 42,31 = 17,31 \quad (\approx 17\,310 \text{ €})
\end{aligned}[/math]

Profitto totale: [math]7\,269[/math] €

Verifica vincoli:
Rischio totale: [math]40(40,38) + 20(42,31) + 60(17,31) = 3500[/math] (Rischio medio 35).
[math]CO_2[/math] risparmiata: [math]500(40,38) + 800(42,31) + 1500(17,31) = 80\,003[/math] kg.

9. Interpretazione

Questa soluzione rappresenta il massimo profitto ottenibile rispettando i vincoli. Si trova sulla superficie di Pareto: ogni miglioramento del profitto richiederebbe di sacrificare il rischio o la sostenibilità.

💡 Osservazione: Il metodo [math]\epsilon[/math]-constraint ha trasformato un problema multi-obiettivo in programmazione lineare classica.

Domanda di riflessione: Immagina di aumentare il requisito di [math]CO_2[/math]. Cosa succederà alla regione ammissibile e al profitto?

Risposta: Aumentando il requisito, la regione si restringe verso il progetto C. Poiché C ha il rendimento più basso, il profitto massimo diminuirà. Questo mostra il costo opportunità della sostenibilità.

Opzione B: Il Metodo “Power” con Python

Useremo Python per generare centinaia di soluzioni ottimali variando sistematicamente i pesi ([math]w_1, w_2, w_3[/math]) assegnati agli obiettivi. Questo ci permetterà di “scansionare” l’intera superficie dei compromessi.

Utilizzeremo la libreria scipy.optimize per trovare i punti di ottimo e matplotlib per la visualizzazione.

Il Codice Python

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# --- 1. Definizione dei Progetti (Dati dell'Esercizio 2) ---
# Budget Totale = 100.000 €
# x[0]=Solare, x[1]=Eolico, x[2]=Idrogeno

def get_performance(x):
    f1_profit = 0.10*x[0] + 0.06*x[1] + 0.04*x[2]    # Profitto (€)
    f2_risk = 40*x[0] + 20*x[1] + 60*x[2]          # Rischio (scala totale)
    f3_co2 = 500*x[0] + 800*x[1] + 1500*x[2]      # CO2 risparmiata (kg)
    return f1_profit, f2_risk, f3_co2

# --- 2. Funzione Obiettivo Scalare (Somma Pesata) ---
# Da minimizzare: invertiamo il segno per Profitto e CO2.
# Normalizziamo Rischio e CO2 per renderli confrontabili con i pesi 0-1.
def scalar_objective(x, w1, w2, w3):
    f1, f2, f3 = get_performance(x)
    return -(w1*(f1/10000) - w2*(f2/6000000) + w3*(f3/150000))

# --- 3. Vincoli e Confini ---
# Vincolo di uguaglianza: Somma investimenti = 100.000
cons = ({'type': 'eq', 'fun': lambda x:  np.sum(x) - 100000})
# Confini: investimenti non negativi
bnds = ((0, None), (0, None), (0, None))

# --- 4. Generazione della Superficie di Pareto ---
pareto_points = []
# Scansioniamo i pesi w1 (Profitto) e w2 (Rischio). w3 (CO2) è il resto.
for w1 in np.linspace(0, 1, 20):
    for w2 in np.linspace(0, 1 - w1, 20):
        w3 = 1 - w1 - w2
        # Risolviamo l'ottimizzazione per questa combinazione di pesi
        res = minimize(scalar_objective, [33333, 33333, 33334], args=(w1, w2, w3), 
                       method='SLSQP', bounds=bnds, constraints=cons)
        if res.success:
            # Salviamo le prestazioni reali del punto trovato
            p, r, c = get_performance(res.x)
            pareto_points.append((p, r, c))

# --- 5. Visualizzazione Tridimensionale su 2D (Colore per CO2) ---
data = np.array(pareto_points)
profit = data[:,0]
risk = data[:,1]
co2 = data[:,2]

plt.figure(figsize=(10, 7))
# Grafico a dispersione: X=Rischio, Y=Profitto, Colore=CO2
scatter = plt.scatter(risk, profit, c=co2, cmap='viridis', alpha=0.7, edgecolors='k')

# Aggiungiamo una barra dei colori per spiegare la CO2
cbar = plt.colorbar(scatter)
cbar.set_label('CO2 risparmiata (kg/anno)', rotation=270, labelpad=15)

plt.title("Mappa dei Compromessi: Profitto vs Rischio vs CO2")
plt.xlabel("Rischio Totale Portafoglio")
plt.ylabel("Profitto Totale Previsionale (€)")
plt.grid(True)
plt.show()