I have the following code that plots the results of two mathematical MILPS with different parameter combinations and determines the Pareto front based on the points. The code is:

import pandas as pd
import matplotlib.pyplot as plt

# Data
data = {
    'understaff1': [0.1, 0.2, 0.15, 0.3, 0.25],
    'viol1': [0.4, 0.35, 0.45, 0.5, 0.55],
    'omega1': [1, 1, 1, 2, 2],
    'gamma1': [0.1, 0.2, 0.1, 0.2, 0.3],
    'understaff2': [0.15, 0.15, 0.2, 0.25, 0.3],
    'viol2': [0.37, 0.4, 0.35, 0.45, 0.5],
    'omega2': [2, 2, 2, 3, 3],
    'gamma2': [0.1, 0.1, 0.2, 0.2, 0.3]
}

df1 = pd.DataFrame({
    'understaff': data['understaff1'],
    'viol': data['viol1'],
    'omega': data['omega1'],
    'gamma': data['gamma1']
})

df2 = pd.DataFrame({
    'understaff': data['understaff2'],
    'viol': data['viol2'],
    'omega': data['omega2'],
    'gamma': data['gamma2']
})

df = pd.concat([df1, df2], ignore_index=True)


# Pareto-Frontier
def pareto_frontier(df):
    pareto_front = []
    for i, row in df.iterrows():
        dominated = False
        for j, other in df.iterrows():
            if (other['understaff'] <= row['understaff'] and other['viol'] <= row['viol']) and (
                    other['understaff'] < row['understaff'] or other['viol'] < row['viol']):
                dominated = True
                break
        if not dominated:
            pareto_front.append(row)
    pareto_front_df = pd.DataFrame(pareto_front)
    pareto_front_df = pareto_front_df.sort_values(by=['understaff'])
    return pareto_front_df


pareto_df = pareto_frontier(df)

# Plot
plt.figure(figsize=(12, 8))
colors = plt.cm.tab20.colors

labels_dict = {}

# First Model
for i, row in df1.iterrows():
    label = f"$omega={row['omega']}, gamma={row['gamma']}$"
    if label not in labels_dict:
        labels_dict[label] = plt.scatter(row['understaff'], row['viol'], color=colors[i % len(colors)], marker='o', s=100, label=label)
    else:
        plt.scatter(row['understaff'], row['viol'], color=colors[i % len(colors)], marker='o', s=100)

# Secon Model
for i, row in df2.iterrows():
    label = f"$omega={row['omega']}, gamma={row['gamma']}$"
    if label not in labels_dict:
        labels_dict[label] = plt.scatter(row['understaff'], row['viol'], color=colors[i % len(colors)], marker='s', s=100, alpha=0.8, label=label)
    else:
        plt.scatter(row['understaff'], row['viol'], color=colors[i % len(colors)], marker='s', s=100, alpha=0.8)

# Pareto-optimal Points
for i, row in pareto_df.iterrows():
    index_in_df1 = df1[(df1['understaff'] == row['understaff']) &
                       (df1['viol'] == row['viol']) &
                       (df1['omega'] == row['omega']) &
                       (df1['gamma'] == row['gamma'])].index
    index_in_df2 = df2[(df2['understaff'] == row['understaff']) &
                       (df2['viol'] == row['viol']) &
                       (df2['omega'] == row['omega']) &
                       (df2['gamma'] == row['gamma'])].index

    if not index_in_df1.empty:
        color_index = index_in_df1[0] % len(colors)
        plt.scatter(row['understaff'], row['viol'], color=colors[color_index], edgecolors='black',
                    linewidths=2, alpha=0.6, s=150, marker='o')
    if not index_in_df2.empty:
        color_index = index_in_df2[0] % len(colors)
        plt.scatter(row['understaff'], row['viol'], color=colors[color_index], edgecolors='black',
                    linewidths=2, alpha=0.6, s=150, marker='s')

# Draw Pareto-Line
plt.plot(pareto_df['understaff'], pareto_df['viol'], linestyle='-', marker='x', color='red', alpha=0.7)

plt.legend(loc='center left', bbox_to_anchor=(1.02, 0.5), ncol=1)

plt.xlabel('Scaled Understaffing')
plt.ylabel('Scaled Violation')
plt.grid(True)
plt.tight_layout()
plt.show()

And this is the result:

Please note that the different lists were chosen because the first lists (with 1 at the end) correspond to the results of model 1 and the lists with a 2 at the end correspond to the results of model two. Now I want to solve the model for different instances per parameter combination, and then display the results per parameter combination as a boxplot. I would then display the data something like this.

'understaff1': [[0.1, 0.2], [0.2, 0.1], [0.15, 0.01] [0.3, 0.2], [0.33, 0.25]],
'viol1': [[0.22, 0.32], [0.12, 0.40], [0.31, 0.15] [0.23, 0.12], [0.34, 0.55]],
'omega1': [[1, 2], [2, 2], [3, 2] [1, 3], [3, 1]],
'gamma1': [[0.22, 0.22], [0.02, 0.04], [0.4, 0.25] [0.31, 0.32], [0.34, 0.25]],
'understaff2': [[0.21, 0.12], [0.32, 0.13], [0.13, 0.4] [0.12, 0.23], [0.04, 0.35]],
'viol2': [[0.32, 0.14], [0.12, 0.12], [0.18, 0.15] [0.25, 0.04], [0.32, 0.15]],
'omega2': [[2, 2], [1, 0], [3, 4], [4, 0], [5, 3]],
'gamma2': [[0.2, 0.12], [0.22, 0.14], [0.1, 0.05] [0.35, 0.02], [0.3, 0.05]]

How do I get these double boxplots to be created and the pareto frontier to be calculated based on the mean values of the respective boxplots of the respective parameter combination? This is what I mean by a double boxplot.