Comparacion final entre metodos

Este notebook compara descenso por gradiente, algoritmo evolutivo, PSO y evolucion diferencial sobre los casos trabajados en la parte de optimizacion numerica.

Se incluyen todas las funciones generalizables en 2D y 3D. En el caso de Goldstein-Price y Six-Hump Camel, el analisis se deja en 2D porque en este proyecto se implementaron en su forma clasica bidimensional.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from funciones_gradientes import (
    griewank_gradient,
    goldstein_price_gradient,
    rastrigin_gradient,
    rosenbrock_gradient,
    run_gradient_descent,
    schwefel_gradient,
    six_hump_camel_gradient,
)
from heuristicos.funciones_heuristicas import (
    run_differential_evolution,
    run_evolutionary_algorithm,
    run_particle_swarm_optimization,
)
from funciones_objetivo import (
    griewank,
    goldstein_price,
    rastrigin,
    rosenbrock,
    schwefel,
    six_hump_camel,
)

COMMON_2D_3D = {
    "Rosenbrock": {
        "objective_function": rosenbrock,
        "gradient_function": rosenbrock_gradient,
        "lower_bounds": -2.048,
        "upper_bounds": 2.048,
        "known_optimum_by_dimension": {
            2: np.array([1.0, 1.0], dtype=float),
            3: np.array([1.0, 1.0, 1.0], dtype=float),
        },
        "gradient_rate": 0.001,
        "gradient_iterations": 5000,
    },
    "Rastrigin": {
        "objective_function": rastrigin,
        "gradient_function": rastrigin_gradient,
        "lower_bounds": -5.12,
        "upper_bounds": 5.12,
        "known_optimum_by_dimension": {
            2: np.zeros(2, dtype=float),
            3: np.zeros(3, dtype=float),
        },
        "gradient_rate": 0.001,
        "gradient_iterations": 4000,
    },
    "Schwefel": {
        "objective_function": schwefel,
        "gradient_function": schwefel_gradient,
        "lower_bounds": -500.0,
        "upper_bounds": 500.0,
        "known_optimum_by_dimension": {
            2: np.array([420.968746, 420.968746], dtype=float),
            3: np.array([420.968746, 420.968746, 420.968746], dtype=float),
        },
        "gradient_rate": 0.01,
        "gradient_iterations": 4000,
    },
    "Griewank": {
        "objective_function": griewank,
        "gradient_function": griewank_gradient,
        "lower_bounds": -600.0,
        "upper_bounds": 600.0,
        "known_optimum_by_dimension": {
            2: np.zeros(2, dtype=float),
            3: np.zeros(3, dtype=float),
        },
        "gradient_rate": 0.01,
        "gradient_iterations": 4000,
    },
}

ONLY_2D = {
    "Goldstein-Price 2D": {
        "objective_function": goldstein_price,
        "gradient_function": goldstein_price_gradient,
        "dimension": 2,
        "lower_bounds": -2.0,
        "upper_bounds": 2.0,
        "known_optimum": np.array([0.0, -1.0], dtype=float),
        "gradient_rate": 0.00001,
        "gradient_iterations": 6000,
    },
    "Six-Hump Camel 2D": {
        "objective_function": six_hump_camel,
        "gradient_function": six_hump_camel_gradient,
        "dimension": 2,
        "lower_bounds": -3.0,
        "upper_bounds": 3.0,
        "known_optimum": np.array([0.089842, -0.712656], dtype=float),
        "gradient_rate": 0.01,
        "gradient_iterations": 5000,
    },
}

def build_cases() -> dict:
    cases = {}

    for function_name, config in COMMON_2D_3D.items():
        for dimension in [2, 3]:
            case_name = f"{function_name} {dimension}D"
            cases[case_name] = {
                "objective_function": config["objective_function"],
                "gradient_function": config["gradient_function"],
                "dimension": dimension,
                "lower_bounds": config["lower_bounds"],
                "upper_bounds": config["upper_bounds"],
                "known_optimum": config["known_optimum_by_dimension"][dimension],
                "gradient_rate": config["gradient_rate"],
                "gradient_iterations": config["gradient_iterations"],
                "gradient_tolerance": 1e-8,
                "ea_population_size": 40 if dimension == 2 else 50,
                "ea_iterations": 120 if dimension == 2 else 150,
                "pso_swarm_size": 40 if dimension == 2 else 50,
                "pso_iterations": 120 if dimension == 2 else 150,
                "de_population_size": 40 if dimension == 2 else 50,
                "de_iterations": 120 if dimension == 2 else 150,
                "seed": 42,
            }

    for case_name, config in ONLY_2D.items():
        cases[case_name] = {
            **config,
            "gradient_tolerance": 1e-8,
            "ea_population_size": 40,
            "ea_iterations": 120,
            "pso_swarm_size": 40,
            "pso_iterations": 120,
            "de_population_size": 40,
            "de_iterations": 120,
            "seed": 42,
        }

    return cases


CASES = build_cases()
list(CASES.keys())

def random_initial_position(dimension: int, lower_bounds: float, upper_bounds: float, seed: int) -> np.ndarray:
    rng = np.random.default_rng(seed)
    return rng.uniform(lower_bounds, upper_bounds, size=dimension)


def summarize_result(case_name: str, method_name: str, result: dict, known_optimum: np.ndarray) -> dict:
    if "final_position" in result:
        best_solution = np.array(result["final_position"], dtype=float)
        best_value = float(result["final_value"])
        iterations = int(result["iterations"])
        evaluations = iterations + 1
    else:
        best_solution = np.array(result["best_solution"], dtype=float)
        best_value = float(result["best_value"])
        iterations = int(result["iterations"])
        evaluations = int(result["evaluations"])

    distance_to_optimum = float(np.linalg.norm(best_solution - known_optimum))

    return {
        "caso": case_name,
        "metodo": method_name,
        "mejor_valor": best_value,
        "iteraciones": iterations,
        "evaluaciones": evaluations,
        "distancia_al_optimo": distance_to_optimum,
        "mejor_solucion": np.array2string(best_solution, precision=4),
    }


def ejecutar_caso(case_name: str, config: dict) -> tuple[pd.DataFrame, dict]:
    objective_function = config["objective_function"]
    gradient_function = config["gradient_function"]
    dimension = config["dimension"]
    lower_bounds = config["lower_bounds"]
    upper_bounds = config["upper_bounds"]
    known_optimum = config["known_optimum"]
    seed = config["seed"]

    initial_position = random_initial_position(dimension, lower_bounds, upper_bounds, seed)

    gradient_result = run_gradient_descent(
        initial_position=initial_position,
        function=objective_function,
        gradient_function=gradient_function,
        rate=config["gradient_rate"],
        max_iterations=config["gradient_iterations"],
        tolerance=config["gradient_tolerance"],
    )

    evolutionary_result = run_evolutionary_algorithm(
        objective_function=objective_function,
        population_size=config["ea_population_size"],
        dimension=dimension,
        lower_bounds=lower_bounds,
        upper_bounds=upper_bounds,
        elitism_fraction=0.2,
        mutation_fraction=0.1,
        max_iterations=config["ea_iterations"],
        seed=seed,
    )

    pso_result = run_particle_swarm_optimization(
        objective_function=objective_function,
        swarm_size=config["pso_swarm_size"],
        dimension=dimension,
        lower_bounds=lower_bounds,
        upper_bounds=upper_bounds,
        inertia_weight=0.7,
        cognitive_weight=1.5,
        social_weight=1.5,
        max_iterations=config["pso_iterations"],
        seed=seed,
    )

    de_result = run_differential_evolution(
        objective_function=objective_function,
        population_size=config["de_population_size"],
        dimension=dimension,
        lower_bounds=lower_bounds,
        upper_bounds=upper_bounds,
        mutation_factor=0.8,
        crossover_rate=0.7,
        max_iterations=config["de_iterations"],
        seed=seed,
    )

    results = {
        "Descenso por gradiente": gradient_result,
        "Algoritmo evolutivo": evolutionary_result,
        "PSO": pso_result,
        "Evolucion diferencial": de_result,
    }

    summary_rows = [
        summarize_result(case_name, method_name, result, known_optimum)
        for method_name, result in results.items()
    ]

    summary_df = pd.DataFrame(summary_rows).sort_values(by=["mejor_valor", "distancia_al_optimo"])
    return summary_df, results

resumenes = []
resultados_por_caso = {}

for case_name, config in CASES.items():
    summary_df, results = ejecutar_caso(case_name, config)
    resumenes.append(summary_df)
    resultados_por_caso[case_name] = results

comparacion_df = pd.concat(resumenes, ignore_index=True)
comparacion_df

for case_name in CASES:
    print(case_name)
    display(
        comparacion_df[comparacion_df["caso"] == case_name]
        .sort_values(by=["mejor_valor", "distancia_al_optimo"])
        .reset_index(drop=True)
    )

tabla_discusion = comparacion_df[["caso", "metodo", "mejor_valor", "evaluaciones", "distancia_al_optimo"]].copy()
tabla_discusion = tabla_discusion.sort_values(by=["caso", "mejor_valor", "evaluaciones"])
tabla_discusion

ranking = (
    comparacion_df.copy()
    .sort_values(by=["caso", "mejor_valor", "distancia_al_optimo", "evaluaciones"])
    .groupby("caso")
    .head(1)
    [["caso", "metodo", "mejor_valor", "evaluaciones", "distancia_al_optimo"]]
    .reset_index(drop=True)
)

ranking

fig, axes = plt.subplots(1, 2, figsize=(16, 5))

promedio_valor = comparacion_df.groupby("metodo")["mejor_valor"].mean().sort_values()
promedio_eval = comparacion_df.groupby("metodo")["evaluaciones"].mean().sort_values()

axes[0].bar(promedio_valor.index, promedio_valor.values, color="steelblue")
axes[0].set_title("Promedio del mejor valor final por metodo")
axes[0].set_ylabel("Valor promedio")
axes[0].tick_params(axis="x", rotation=20)
axes[0].grid(axis="y", alpha=0.3)

axes[1].bar(promedio_eval.index, promedio_eval.values, color="indianred")
axes[1].set_title("Promedio de evaluaciones por metodo")
axes[1].set_ylabel("Evaluaciones promedio")
axes[1].tick_params(axis="x", rotation=20)
axes[1].grid(axis="y", alpha=0.3)

plt.tight_layout()
plt.show()

selected_cases = ["Rosenbrock 2D", "Rastrigin 2D", "Schwefel 3D", "Griewank 3D"]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
axes = axes.ravel()

for ax, case_name in zip(axes, selected_cases):
    results = resultados_por_caso[case_name]

    ax.plot(results["Descenso por gradiente"]["function_values"], label="Gradiente", linewidth=2)
    ax.plot(results["Algoritmo evolutivo"]["best_values_history"], label="Alg. evolutivo", linewidth=2)
    ax.plot(results["PSO"]["best_values_history"], label="PSO", linewidth=2)
    ax.plot(results["Evolucion diferencial"]["best_values_history"], label="Evol. diferencial", linewidth=2)
    ax.set_title(f"Convergencia - {case_name}")
    ax.set_xlabel("Iteracion")
    ax.set_ylabel("Mejor valor encontrado")
    ax.set_yscale("symlog")
    ax.grid(alpha=0.3)

axes[0].legend()
plt.tight_layout()
plt.show()

Lectura sugerida de resultados

• mejor_valor: que tan bien minimizo cada metodo.
• evaluaciones: cuantas evaluaciones de la funcion objetivo uso cada metodo.
• distancia_al_optimo: que tan cerca quedo de la solucion conocida.
• ranking: resume el mejor metodo por caso con base en valor final, cercania al optimo y evaluaciones.
• Las curvas de convergencia se muestran solo para algunos casos representativos para no saturar la visualizacion final.

Nota para el reporte

En la redaccion final conviene aclarar que Goldstein-Price y Six-Hump Camel se dejaron en 2D porque en el proyecto se implementaron en su forma clasica bidimensional. Asi la comparacion queda completa para todos los casos efectivamente desarrollados y no genera una inconsistencia metodologica con el codigo.