---

from dataclasses import dataclass
import math
import random
import time

import numpy as np
import matplotlib.pyplot as plt


# ============================================================
# 1. Параметры задачи и ГА
# ============================================================

SEED = 42
POPULATION_SIZE = 80
GENERATIONS = 200
ELITISM = 4
TOURNAMENT_SIZE = 3
CROSSOVER_RATE = 0.90
MUTATION_RATE = 0.10
STOP_WITHOUT_IMPROVEMENT = 60


E = 2.1e11                # Па
DENSITY = 7850.0          # кг/м³
SIGMA_ALLOW = 160e6       # Па
DISPLACEMENT_ALLOW = 0.020  # м (20 мм)

# Геометрические ограничения
MIN_MEMBER_LENGTH = 0.10
MIN_NODE_DISTANCE = 0.10
MU_BUCKLING = 1.0
BUCKLING_SAFETY_FACTOR = 1.0

# Сортамент: КРУГЛЫЕ СТЕРЖНИ СПЛОШНОГО СЕЧЕНИЯ, диаметр 12–100 мм.
# A = πd²/4, I = πd⁴/64 — формулы сплошного круга.
DIAMETERS_MM = np.array([
    12, 14, 16, 18, 20, 22, 25, 28, 30, 32,
    35, 38, 40, 45, 50, 55, 60, 70, 80, 90, 100
], dtype=float)

AREAS = np.pi * (DIAMETERS_MM / 1000.0) ** 2 / 4.0
INERTIAS = np.pi * (DIAMETERS_MM / 1000.0) ** 4 / 64.0


# ============================================================
# 2. Тестовая ферма
# ============================================================

NODES = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [2.0, 0.0],
    [0.0, 1.0],
    [1.0, 1.0],
    [2.0, 1.0],
], dtype=float)

MEMBERS = np.array([
    [0, 1], [1, 2],
    [3, 4], [4, 5],
    [0, 3], [1, 4], [2, 5],
    [0, 4], [1, 3], [1, 5],
], dtype=int)

SUPPORTS = {
    0: (True, True),
    2: (False, True),
}

LOADS = np.zeros((len(NODES), 2), dtype=float)
LOADS[1, 1] = -30e3
LOADS[4, 1] = -30e3

SHAPE_BOUNDS = {
    3: ((0.0, 0.0), (0.6, 1.5)),
    4: ((0.6, 1.4), (0.6, 1.5)),
    5: ((2.0, 2.0), (0.6, 1.5)),
}


# ============================================================
# 3. Структура результата
# ============================================================

@dataclass
class Result:
    fitness: float
    mass: float
    feasible: bool
    violations: float
    coords: np.ndarray
    diameters_mm: np.ndarray
    active: np.ndarray
    displacements: np.ndarray
    axial_forces: np.ndarray
    stresses: np.ndarray
    buckling_ratios: np.ndarray
    reason: str = ""


# ============================================================
# 4. МКЭ
# ============================================================

def verify_rigidity(coords=None):
    """Проверка геометрической неизменяемости исходной схемы."""
    c = NODES if coords is None else coords
    n = len(c)
    ndof = 2 * n
    K = np.zeros((ndof, ndof))
    areas = np.ones(len(MEMBERS)) * 1e-4
    for e, (i, j) in enumerate(MEMBERS):
        d = c[j] - c[i]
        L = float(np.linalg.norm(d))
        if L < 1e-9:
            raise ValueError(f"Стержень {e} нулевой длины.")
        cx, cy = d / L
        k = (E * areas[e] / L) * np.array([
            [ cx*cx,  cx*cy, -cx*cx, -cx*cy],
            [ cx*cy,  cy*cy, -cx*cy, -cy*cy],
            [-cx*cx, -cx*cy,  cx*cx,  cx*cy],
            [-cx*cy, -cy*cy,  cx*cy,  cy*cy],
        ])
        dofs = [2*i, 2*i+1, 2*j, 2*j+1]
        K[np.ix_(dofs, dofs)] += k

    fixed = set()
    for node, (fx, fy) in SUPPORTS.items():
        if fx: fixed.add(2 * node)
        if fy: fixed.add(2 * node + 1)
    free = np.array([d for d in range(ndof) if d not in fixed], dtype=int)

    K_ff = K[np.ix_(free, free)]
    rank = np.linalg.matrix_rank(K_ff)
    cond = np.linalg.cond(K_ff)
    return rank, K_ff.shape[0], cond


def assemble_and_solve(coords, areas, active):

    n = len(coords)
    ndof = 2 * n
    K = np.zeros((ndof, ndof), dtype=float)
    lengths = np.zeros(len(MEMBERS), dtype=float)
    dirs = np.zeros((len(MEMBERS), 2), dtype=float)

    for e, (i, j) in enumerate(MEMBERS):
        if not active[e]:
            continue
        d = coords[j] - coords[i]
        L = float(np.linalg.norm(d))
        if L < 1e-9:
            raise ValueError(f"Стержень {e} нулевой длины.")
        cx, cy = d / L
        lengths[e] = L
        dirs[e] = (cx, cy)
        k = (E * areas[e] / L) * np.array([
            [ cx*cx,  cx*cy, -cx*cx, -cx*cy],
            [ cx*cy,  cy*cy, -cx*cy, -cy*cy],
            [-cx*cx, -cx*cy,  cx*cx,  cx*cy],
            [-cx*cy, -cy*cy,  cx*cy,  cy*cy],
        ])
        dofs = [2*i, 2*i+1, 2*j, 2*j+1]
        K[np.ix_(dofs, dofs)] += k

    F = LOADS.reshape(-1)
    fixed = set()
    for node, (fx, fy) in SUPPORTS.items():
        if fx: fixed.add(2 * node)
        if fy: fixed.add(2 * node + 1)
    free = np.array([d for d in range(ndof) if d not in fixed], dtype=int)
    if len(free) == 0:
        raise ValueError("Нет свободных DOF.")

    K_ff = K[np.ix_(free, free)]
    F_f = F[free]

    try:
        u_free = np.linalg.solve(K_ff, F_f)
    except np.linalg.LinAlgError as exc:
        raise ValueError("Сингулярная матрица.") from exc

    if not np.all(np.isfinite(u_free)):
        raise ValueError("NaN/Inf в решении.")
    # Защита от катастрофически больших чисел (почти-сингулярная матрица)
    if np.max(np.abs(u_free)) > 1e3:
        raise ValueError("Решение аномально большое (почти-сингулярная матрица).")

    u = np.zeros(ndof)
    u[free] = u_free

    forces = np.zeros(len(MEMBERS))
    stresses = np.zeros(len(MEMBERS))
    for e, (i, j) in enumerate(MEMBERS):
        if not active[e]:
            continue
        L = lengths[e]
        cx, cy = dirs[e]
        ext = np.dot(np.array([-cx, -cy, cx, cy]),
                     u[[2*i, 2*i+1, 2*j, 2*j+1]])
        forces[e] = E * areas[e] * ext / L
        stresses[e] = forces[e] / areas[e]

    return u.reshape((-1, 2)), forces, stresses, lengths


# ============================================================
# 5. Кодирование и оценка
# ============================================================

def random_chromosome():
    genes = []
    for _ in MEMBERS:
        genes.append(random.randrange(len(DIAMETERS_MM)))
    for node, bounds in SHAPE_BOUNDS.items():
        for low, high in bounds:
            genes.append(random.uniform(low, high))
    return genes


def initial_chromosome():
    genes = []
    idx_25 = int(np.argmin(np.abs(DIAMETERS_MM - 25.0)))
    for _ in MEMBERS:
        genes.append(idx_25)
    for node, bounds in SHAPE_BOUNDS.items():
        base = NODES[node]
        for axis, (low, high) in enumerate(bounds):
            genes.append(float(np.clip(base[axis], low, high)))
    return genes


def decode(chromosome):
    n = len(MEMBERS)
    idx = np.asarray(chromosome[:n], dtype=int)
    diameters = DIAMETERS_MM[np.clip(idx, 0, len(DIAMETERS_MM) - 1)]
    coords = NODES.copy()
    offset = n
    for node, bounds in SHAPE_BOUNDS.items():
        for axis, (low, high) in enumerate(bounds):
            coords[node, axis] = float(np.clip(chromosome[offset], low, high))
            offset += 1
    active = np.ones(n, dtype=bool)
    return coords, diameters, active


def evaluate(chromosome):
    coords, diameters, active = decode(chromosome)
    areas = np.zeros(len(MEMBERS))
    inertias = np.zeros(len(MEMBERS))
    for e, d in enumerate(diameters):
        k = int(np.argmin(np.abs(DIAMETERS_MM - d)))
        areas[e] = AREAS[k]
        inertias[e] = INERTIAS[k]

    # Геометрические проверки
    for e, (i, j) in enumerate(MEMBERS):
        L = float(np.linalg.norm(coords[j] - coords[i]))
        if L < MIN_MEMBER_LENGTH:
            return _invalid(coords, diameters, active,
                            f"Стержень {e} короче {MIN_MEMBER_LENGTH} м.")
    for i in range(len(coords)):
        for j in range(i + 1, len(coords)):
            d = float(np.linalg.norm(coords[j] - coords[i]))
            if d < MIN_NODE_DISTANCE:
                return _invalid(coords, diameters, active,
                                f"Узлы {i} и {j} слишком близко.")

    # МКЭ
    try:
        u, forces, stresses, lengths = assemble_and_solve(coords, areas, active)
    except (ValueError, FloatingPointError) as exc:
        return _invalid(coords, diameters, active, str(exc))

    # Масса
    mass = float(np.sum(DENSITY * areas[active] * lengths[active]))

    ratios = []

    # 1) Прочность
    for e in range(len(MEMBERS)):
        if not active[e]:
            continue
        r = max(0.0, abs(stresses[e]) / SIGMA_ALLOW - 1.0)
        ratios.append(r)

    # 2) Результирующие перемещения sqrt(ux² + uy²) для свободных узлов
    for node in range(len(NODES)):
        free_x = node not in SUPPORTS or not SUPPORTS[node][0]
        free_y = node not in SUPPORTS or not SUPPORTS[node][1]
        if free_x or free_y:
            mag = math.hypot(u[node, 0], u[node, 1])
            r = max(0.0, mag / DISPLACEMENT_ALLOW - 1.0)
            ratios.append(r)

    # 3) Устойчивость сжатых стержней
    buckling = np.zeros(len(MEMBERS))
    for e in range(len(MEMBERS)):
        if not active[e] or forces[e] >= 0:
            continue
        L_eff = MU_BUCKLING * lengths[e]
        p_cr = math.pi**2 * E * inertias[e] / max(L_eff**2, 1e-20)
        p_dem = abs(forces[e]) * BUCKLING_SAFETY_FACTOR
        buckling[e] = p_dem / max(p_cr, 1e-20)
        ratios.append(max(0.0, buckling[e] - 1.0))

    # 4) Эвристика: штраф за недогрузку.
    #    НЕ входит в violation (не влияет на допустимость),
    #    только мягко увеличивает fitness.
    underload = 0.0
    for e in range(len(MEMBERS)):
        if not active[e]:
            continue
        N_max = areas[e] * SIGMA_ALLOW
        util = abs(forces[e]) / max(N_max, 1e-12)
        if util < 0.05:
            underload += (0.05 - util) / 0.05

    violation = float(sum(ratios))
    feasible = violation <= 1e-8

    fitness = (mass
               * (1.0 + 1e4 * violation + 1e4 * violation ** 2)
               * (1.0 + 0.2 * underload))

    return Result(
        fitness=fitness, mass=mass, feasible=feasible,
        violations=violation, coords=coords,
        diameters_mm=diameters, active=active,
        displacements=u, axial_forces=forces,
        stresses=stresses, buckling_ratios=buckling, reason=""
    )


def _invalid(coords, diameters, active, reason):
    return Result(
        fitness=1e12, mass=1e9, feasible=False, violations=1e6,
        coords=coords, diameters_mm=diameters, active=active,
        displacements=np.zeros_like(coords),
        axial_forces=np.zeros(len(MEMBERS)),
        stresses=np.zeros(len(MEMBERS)),
        buckling_ratios=np.zeros(len(MEMBERS)),
        reason=reason
    )


# ============================================================
# 6. Генетические операторы
# ============================================================

def tournament(population, scores):
    k = min(TOURNAMENT_SIZE, len(population))
    idx = random.sample(range(len(population)), k)
    return population[min(idx, key=lambda i: scores[i])][:]


def crossover(p1, p2):
    if random.random() > CROSSOVER_RATE or len(p1) < 2:
        return p1[:], p2[:]
    a, b = p1[:], p2[:]
    l, r = sorted(random.sample(range(len(a) + 1), 2))
    a[l:r], b[l:r] = b[l:r], a[l:r]
    return a, b


def mutate(chrom, gen):
    child = chrom[:]
    n = len(MEMBERS)
    progress = gen / max(GENERATIONS - 1, 1)
    step = 0.20 * (1.0 - progress) + 0.02
    for g in range(len(child)):
        if random.random() >= MUTATION_RATE:
            continue
        if g < n:
            child[g] = random.randrange(len(DIAMETERS_MM))
        else:
            si = g - n
            node_pos = si // 2
            axis = si % 2
            node = list(SHAPE_BOUNDS.keys())[node_pos]
            low, high = SHAPE_BOUNDS[node][axis]
            if high > low:
                child[g] = float(np.clip(
                    child[g] + random.gauss(0, step * (high - low)),
                    low, high))
    return child


def invert(chrom):
    child = chrom[:]
    if len(child) >= 4 and random.random() < 0.10:
        l, r = sorted(random.sample(range(len(child)), 2))
        child[l:r + 1] = child[l:r + 1][::-1]
    return child


# ============================================================
# 7. Основной цикл
# ============================================================

def optimize(verbose=True):
    random.seed(SEED)
    np.random.seed(SEED)

    population = [initial_chromosome()] + \
                 [random_chromosome() for _ in range(POPULATION_SIZE - 1)]

    best_feasible_chrom = None
    best_feasible_mass = float("inf")
    best_seen = None
    best_score = float("inf")
    stagnant = 0
    history_best = []
    history_mean = []
    history_best_mass = []
    history_feasible_frac = []
    t0 = time.time()

    for gen in range(GENERATIONS):
        results = [evaluate(ind) for ind in population]
        scores = [r.fitness for r in results]

        order = sorted(range(len(population)), key=lambda i: scores[i])
        population = [population[i] for i in order]
        results = [results[i] for i in order]
        scores = [scores[i] for i in order]

        history_best.append(scores[0])
        history_mean.append(float(np.mean(scores)))
        feasible_frac = float(np.mean([r.feasible for r in results]))
        history_feasible_frac.append(feasible_frac)

        for ind, r in zip(population, results):
            if r.feasible and r.mass < best_feasible_mass:
                best_feasible_mass = r.mass
                best_feasible_chrom = ind[:]

        history_best_mass.append(best_feasible_mass
                                 if best_feasible_chrom is not None
                                 else float("nan"))

        current = results[0]
        if current.fitness < best_score:
            best_score = current.fitness
            best_seen = population[0][:]
            stagnant = 0
        else:
            stagnant += 1

        if verbose and (gen % 20 == 0 or gen == GENERATIONS - 1):
            if best_feasible_chrom is not None:
                print(f"Пок.{gen:4d} | допустимых "
                      f"{int(feasible_frac*POPULATION_SIZE):2d}/{POPULATION_SIZE} "
                      f"| лучшая допустимая масса = {best_feasible_mass:7.3f} кг")
            else:
                print(f"Пок.{gen:4d} | допустимых 0/{POPULATION_SIZE} "
                      f"| f* = {current.fitness:.3e} | "
                      f"наруш. = {current.violations:.3e}")

        if stagnant >= STOP_WITHOUT_IMPROVEMENT:
            if verbose:
                print(f"Остановка: нет улучшения {stagnant} поколений.")
            break

        next_pop = population[:ELITISM]
        while len(next_pop) < POPULATION_SIZE:
            p1 = tournament(population, scores)
            p2 = tournament(population, scores)
            c1, c2 = crossover(p1, p2)
            c1 = invert(mutate(c1, gen))
            c2 = invert(mutate(c2, gen))
            next_pop.append(c1)
            if len(next_pop) < POPULATION_SIZE:
                next_pop.append(c2)
        population = next_pop

    elapsed = time.time() - t0

    if best_feasible_chrom is not None:
        best = evaluate(best_feasible_chrom)
    else:
        best = evaluate(best_seen)

    return best, elapsed, history_best, history_mean, history_best_mass, history_feasible_frac


# ============================================================
# 8. Визуализация
# ============================================================

def plot_truss(coords, members, forces=None, ax=None, title="Ферма"):
    if ax is None:
        _, ax = plt.subplots(figsize=(8, 5))
    for e, (i, j) in enumerate(members):
        x = [coords[i, 0], coords[j, 0]]
        y = [coords[i, 1], coords[j, 1]]
        if forces is not None and np.isfinite(forces[e]):
            color = 'red' if forces[e] > 0 else 'blue'
            lw = 1.0 + 3.0 * min(abs(forces[e]) / 1e5, 1.0)
        else:
            color = 'gray'; lw = 1.5
        ax.plot(x, y, color=color, lw=lw, zorder=1)
    ax.scatter(coords[:, 0], coords[:, 1], c='black', s=30, zorder=3)
    for i, (x, y) in enumerate(coords):
        ax.annotate(str(i), (x, y), textcoords="offset points",
                    xytext=(5, 5), fontsize=9)
    for node, (fx, fy) in SUPPORTS.items():
        x, y = coords[node]
        marker = '^' if (fx and fy) else 'o'
        ax.plot(x, y, marker=marker, markersize=12, color='green', zorder=4)
    for node in range(len(coords)):
        Fx, Fy = LOADS[node]
        if Fx == 0 and Fy == 0:
            continue
        x, y = coords[node]
        s = 5e-5
        ax.arrow(x, y, Fx*s, Fy*s, head_width=0.05,
                 head_length=0.05, fc='orange', ec='orange', zorder=5)
    ax.set_aspect('equal')
    ax.set_title(title)
    ax.grid(True, alpha=0.3)


def plot_convergence(hist_best, hist_mean, hist_mass, hist_frac):
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    ax = axes[0]
    ax.plot(hist_best, label='fitness (лучшее)', lw=2)
    ax.set_yscale('log')
    ax.set_xlabel('Поколение')
    ax.set_ylabel('Fitness')
    ax.set_title('Сходимость fitness')
    ax.legend()
    ax.grid(True, alpha=0.3)

    ax = axes[1]
    ax.plot(hist_mass, label='Лучшая допустимая масса, кг', lw=2, color='green')
    ax.set_xlabel('Поколение')
    ax.set_ylabel('Масса, кг')
    ax.set_title('Лучшая допустимая масса')
    ax2 = ax.twinx()
    ax2.plot(hist_frac, label='Доля допустимых', color='orange',
             lw=1.2, alpha=0.7)
    ax2.set_ylabel('Доля допустимых', color='orange')
    ax.legend(loc='upper right')
    ax2.legend(loc='lower right')
    ax.grid(True, alpha=0.3)

    plt.tight_layout()
    return fig


# ============================================================
# 9. Вывод результата
# ============================================================

def print_result(r, elapsed, mass0=None):
    print("\n" + "=" * 72)
    print(" РЕЗУЛЬТАТ ОПТИМИЗАЦИИ")
    print("=" * 72)
    if mass0 is not None:
        print(f"  Масса исходная:       {mass0:.3f} кг")
    print(f"  Масса оптимальная:    {r.mass:.3f} кг")
    if mass0 is not None and mass0 > 0:
        print(f"  Снижение массы:       {(1 - r.mass/mass0)*100:.1f}%")
    print(f"  Допустимое решение:   {'ДА' if r.feasible else 'НЕТ'}")
    print(f"  Нарушения:            {r.violations:.4e}")
    print(f"  Время:                {elapsed:.2f} с")

    print(f"\n  Координаты узлов, м:")
    for i, xy in enumerate(r.coords):
        print(f"    Узел {i}: x = {xy[0]:.4f}, y = {xy[1]:.4f}")

    print(f"\n  Стержни:")
    print(f"    {'№':>3s}  {'(i,j)':>8s}  {'d, мм':>7s}  "
          f"{'N, кН':>9s}  {'σ, МПа':>9s}  {'Эйлер':>8s}")
    for e, (i, j) in enumerate(MEMBERS):
        print(f"    {e:>3d}  {f'({i},{j})':>8s}  {r.diameters_mm[e]:7.2f}  "
              f"{r.axial_forces[e]/1e3:+9.3f}  "
              f"{r.stresses[e]/1e6:+9.2f}  "
              f"{r.buckling_ratios[e]:8.4f}")

    print(f"\n  Перемещения узлов, мм (результирующее):")
    for i, uv in enumerate(r.displacements):
        mag = math.hypot(uv[0], uv[1]) * 1000
        print(f"    Узел {i}: ux = {uv[0]*1000:+8.3f}, "
              f"uy = {uv[1]*1000:+8.3f}, |u| = {mag:7.3f}")

    print(f"\n  Проверка ограничений:")
    print(f"    Max |σ| / σ_доп = {np.max(np.abs(r.stresses))/SIGMA_ALLOW:.3f}")
    u_mag = np.max(np.hypot(r.displacements[:, 0], r.displacements[:, 1]))
    print(f"    Max |u| / u_доп = {u_mag/DISPLACEMENT_ALLOW:.3f}")
    print(f"    Max Эйлер       = {np.max(r.buckling_ratios):.3f}")

    if not r.feasible:
        print(f"\n  ВНИМАНИЕ: решение недопустимо. {r.reason}")


# ============================================================
# 10. Главная функция
# ============================================================

def main():
    print("=" * 72)
    print(" ОПТИМИЗАЦИЯ ПЛОСКОЙ ФЕРМЫ: МКЭ + ГА  (v16)")
    print("=" * 72)

    # Проверка жёсткости исходной схемы
    print("\n  Проверка геометрической неизменяемости исходной схемы:")
    try:
        rank, n_dof, cond = verify_rigidity(NODES)
        print(f"    Свободных DOF: {n_dof}, ранг K_ff: {rank}, "
              f"cond = {cond:.2e}")
        if rank < n_dof:
            print("    !!! СХЕМА ГЕОМЕТРИЧЕСКИ ИЗМЕНЯЕМА !!!")
            return
        print("    Схема геометрически неизменяема.")
    except ValueError as exc:
        print(f"    Ошибка проверки: {exc}")
        return

    print(f"\n  Материал: сталь С245")
    print(f"  E = {E/1e9:.0f} ГПа, ρ = {DENSITY} кг/м³, "
          f"σ_доп = {SIGMA_ALLOW/1e6:.0f} МПа")
    print(f"  Сортамент: круглые стержни сплошного сечения, "
          f"d = {DIAMETERS_MM[0]:.0f}–{DIAMETERS_MM[-1]:.0f} мм")
    print(f"  Узлов: {len(NODES)}, стержней: {len(MEMBERS)}, "
          f"опорных связей: {sum(sum(v) for v in SUPPORTS.values())}")
    m = len(MEMBERS); r_ = sum(sum(v) for v in SUPPORTS.values())
    j = len(NODES)
    print(f"  W = m + r - 2j = {m} + {r_} - {2*j} = {m + r_ - 2*j} "
          f"(степень статической неопределимости)")

    # Исходная ферма: d = 25 мм
    d0 = 25.0
    k0 = int(np.argmin(np.abs(DIAMETERS_MM - d0)))
    areas0 = AREAS[k0] * np.ones(len(MEMBERS))
    active0 = np.ones(len(MEMBERS), dtype=bool)
    try:
        u0, f0, s0, L0 = assemble_and_solve(NODES, areas0, active0)
        mass0 = float(np.sum(DENSITY * areas0 * L0))
        print(f"\n  Исходная ферма (все d = {d0:.0f} мм):")
        print(f"    Масса: {mass0:.3f} кг")
        u0_mag = np.max(np.hypot(u0[:, 0], u0[:, 1]))
        print(f"    Max |u| = {u0_mag*1000:.3f} мм")
        print(f"    Max |σ| = {np.max(np.abs(s0))/1e6:.2f} МПа")
        ratios = []
        for e in range(len(MEMBERS)):
            if f0[e] < 0:
                L_eff = MU_BUCKLING * L0[e]
                p_cr = math.pi**2 * E * INERTIAS[k0] / L_eff**2
                ratios.append(abs(f0[e]) / p_cr)
        print(f"    Max Эйлер = {max(ratios) if ratios else 0:.3f}")
    except ValueError as exc:
        print(f"\n  Исходная ферма не решается: {exc}")
        return

    print("\n  Запуск генетического алгоритма...")
    (result, elapsed, hb, hm, hbm, hbf) = optimize(verbose=True)
    print_result(result, elapsed, mass0)

    # Визуализация
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    plot_truss(NODES, MEMBERS, forces=f0,
               ax=axes[0], title=f"Исходная ферма (M = {mass0:.2f} кг)")
    plot_truss(result.coords, MEMBERS, forces=result.axial_forces,
               ax=axes[1],
               title=f"Оптимальная ферма (M = {result.mass:.2f} кг)")
    plt.tight_layout()
    plt.savefig('truss_result.png', dpi=120)

    plot_convergence(hb, hm, hbm, hbf)
    plt.savefig('convergence.png', dpi=120)
    print("\n  Графики: truss_result.png, convergence.png")
    plt.show()


if __name__ == "__main__":
    main()