From e14869691b3250ab1e939a5728d7f34ac0191928 Mon Sep 17 00:00:00 2001 From: Nico Date: Fri, 6 Feb 2026 11:24:13 +0100 Subject: [PATCH] =?UTF-8?q?turbom=C3=A4=C3=9Fige=20Anpassungen?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GHA_triaxial/gha2_ES.py | 48 ++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/GHA_triaxial/gha2_ES.py b/GHA_triaxial/gha2_ES.py index 3ed29aa..b1e778f 100644 --- a/GHA_triaxial/gha2_ES.py +++ b/GHA_triaxial/gha2_ES.py @@ -2,11 +2,12 @@ import numpy as np from Hansen_ES_CMA import escma from ellipsoide import EllipsoidTriaxial from numpy.typing import NDArray +from typing import Tuple import plotly.graph_objects as go from GHA_triaxial.gha2_num import gha2_num -from GHA_triaxial.utils import sigma2alpha +from GHA_triaxial.utils import sigma2alpha, pq_ell + -ell_ES: EllipsoidTriaxial = None P_left: NDArray = None P_right: NDArray = None @@ -20,6 +21,7 @@ def Sehne(P1: NDArray, P2: NDArray) -> float: """ R12 = P2-P1 s = np.linalg.norm(R12) + return s @@ -31,45 +33,42 @@ def midpoint_fitness(x: tuple) -> float: :param x: enthält die Startwerte von u und v :return: Fitnesswert (f) """ - global ell_ES, P_left, P_right + global P_left, P_right u, v = x - P_middle = ell_ES.para2cart(u, v) + P_middle = ell.para2cart(u, v) d1 = Sehne(P_left, P_middle) d2 = Sehne(P_middle, P_right) base = d1 + d2 # midpoint penalty (dimensionslos) - # relative Differenz, skaliert stabil über verschiedene Segmentlängen + # relative Differenz, skaliert über verschiedene Segmentlängen denom = max(base, 1e-9) pen_equal = ((d1 - d2) / denom) ** 2 w_equal = 10.0 f = base + denom * w_equal * pen_equal + return f -def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, maxIter: int = 10000, all_points: bool = False): +def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, all_points: bool = True) -> Tuple[float, float, float, NDArray]: """ Berechnen der 2. GHA mithilfe der CMA-ES. Die CMA-ES optimiert sukzessive den Mittelpunkt zwischen Start- und Zielpunkt. Der Abbruch der Berechnung erfolgt, wenn alle Segmentlängen <= maxSegLen sind. Die Distanzen zwischen den einzelnen Punkten werden als direkte 3D-Distanzen berechnet und aufaddiert. - :param ell: Parameter des triaxialen Ellipsoids + :param ell: Ellipsoid :param P0: Startpunkt :param Pk: Zielpunkt :param maxSegLen: maximale Segmentlänge - :param maxIter: maximale Durchläufe der Mittelpunktsgenerierung - :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben werden kann - :return: Richtungswinkel des Start- und Zielpunktes und Gesamtlänge + :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben wird + :return: Richtungswinkel in RAD des Start- und Zielpunktes und Gesamtlänge """ - global ell_ES - ell_ES = ell R0 = (ell.ax + ell.ay + ell.b) / 3 if maxSegLen is None: maxSegLen = R0 * 1 / (637.4*2) # 10km Segment bei mittleren Erdradius sigma_uv_nom = 1e-3 * (maxSegLen / R0) # ~1e-5 - points: list[NDArray] = [P0, Pk] startIter = 0 level = 0 @@ -82,7 +81,6 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = level += 1 new_points: list[NDArray] = [points[0]] - for i in range(len(points) - 1): A = points[i] B = points[i+1] @@ -92,22 +90,22 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = if dAB > maxSegLen: global P_left, P_right P_left, P_right = A, B - Au, Av = ell_ES.cart2para(A) - Bu, Bv = ell_ES.cart2para(B) + Au, Av = ell.cart2para(A) + Bu, Bv = ell.cart2para(B) u0 = (Au + Bu) / 2 v0 = Av + 0.5 * np.arctan2(np.sin(Bv - Av), np.cos(Bv - Av)) xmean = [u0, v0] sigmaStep = sigma_uv_nom * (Sehne(A, B) / maxSegLen) - u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep) + u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep) # Aufruf CMA-ES P_next = ell.para2cart(u, v) new_points.append(P_next) startIter += 1 + maxIter = 10000 if startIter > maxIter: raise RuntimeError("Abbruch: maximale Iterationen überschritten.") - new_points.append(B) points = new_points @@ -117,13 +115,13 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = totalLen = float(np.sum(np.linalg.norm(P_all[1:] - P_all[:-1], axis=1))) if len(points) >= 3: - p0i = ell_ES.point_onto_ellipsoid(P0 + 10.0 * (points[1] - P0) / np.linalg.norm(points[1] - P0)) + p0i = ell.point_onto_ellipsoid(P0 + 10.0 * (points[1] - P0) / np.linalg.norm(points[1] - P0)) sigma0 = (p0i - P0) / np.linalg.norm(p0i - P0) - alpha0 = sigma2alpha(ell_ES, sigma0, P0) + alpha0 = sigma2alpha(ell, sigma0, P0) - p1i = ell_ES.point_onto_ellipsoid(Pk - 10.0 * (Pk - points[-2]) / np.linalg.norm(Pk - points[-2])) + p1i = ell.point_onto_ellipsoid(Pk - 10.0 * (Pk - points[-2]) / np.linalg.norm(Pk - points[-2])) sigma1 = (Pk - p1i) / np.linalg.norm(Pk - p1i) - alpha1 = sigma2alpha(ell_ES, sigma1, Pk) + alpha1 = sigma2alpha(ell, sigma1, Pk) else: alpha0 = None alpha1 = None @@ -166,7 +164,7 @@ if __name__ == '__main__': beta0, lamb0 = (0.2, 0.1) P0 = ell.ell2cart(beta0, lamb0) - beta1, lamb1 = (0.7, 0.3) + beta1, lamb1 = (0.3, 0.2) P1 = ell.ell2cart(beta1, lamb1) alpha0, alpha1, s_num, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000, all_points=True) @@ -175,8 +173,10 @@ if __name__ == '__main__': points_num.append(ell.ell2cart(beta, lamb)) points_num = np.array(points_num) - alpha0, alpha1, s, points = gha2_ES(ell, P0, P1, all_points=True) + alpha0, alpha1, s, points = gha2_ES(ell, P0, P1) print(s_num) print(s) + print(alpha0) + print(alpha1) print(s - s_num) show_points(points, points_num, P0, P1)