Pythonfiles

2025-12-28 17:47:23 +01:00
parent e21d39844b
commit cd0d60093e
4 changed files with 239 additions and 107 deletions
--- a/Campusnetz.ipynb
+++ b/Campusnetz.ipynb
@@ -19,7 +19,8 @@
    "import Stochastisches_Modell\n",
    "from Stochastisches_Modell import StochastischesModell\n",
    "import Export\n",
-    "import Netzqualität_Genauigkeit"
+    "import Netzqualität_Genauigkeit\n",
+    "import Datumsfestlegung"
   ],
   "outputs": [],
   "execution_count": null
@@ -418,6 +419,24 @@
   "id": "3b38998a2765b74d",
   "outputs": [],
   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "outputs": [],
+   "execution_count": null,
+   "source": [
+    "# Datumsfestlegung: Bitte geben Sie nachfolgend die Koordinatenkomponenten an, die das Datum definieren sollen\n",
+    "\n",
+    "auswahl = [\n",
+    "    (\"101\",\"X\"), (\"101\",\"Y\"),          # Punkt 101 nur Lage\n",
+    "    (\"205\",\"X\"), (\"205\",\"Y\"), (\"205\",\"Z\"),  # Punkt 205 voll\n",
+    "    (\"330\",\"Z\")                        # Punkt 330 nur Höhe\n",
+    "]\n",
+    "\n",
+    "aktive_unbekannte_indices = Datumsfestlegung.datumskomponenten(auswahl, liste_punktnummern)"
+   ],
+   "id": "c37670a07848d977"
  }
 ],
 "metadata": {
--- a/Datumsfestlegung.py
+++ b/Datumsfestlegung.py
@@ -1,81 +1,154 @@
 import sympy as sp
-from typing import List, Iterable, Tuple
+from typing import Iterable, List, Sequence, Tuple, Optional

-def raenderungsmatrix_G(
-    x0: sp.Matrix,
-    idx_X: List[int],
-    idx_Y: List[int],
-    idx_Z: List[int],
-    mit_massstab: bool = True,
+
+class Datumsfestlegung:
+
+    @staticmethod
+    def datumskomponenten(
+        auswahl: Iterable[Tuple[str, str]],
+        liste_punktnummern: Sequence[str],
+        *,
+        layout: str = "XYZ"
+    ) -> List[int]:
+        punkt2pos = {str(p): i for i, p in enumerate(liste_punktnummern)}
+
+        layout = layout.upper()
+        if layout != "XYZ":
+            raise ValueError("Nur layout='XYZ' unterstützt (wie bei euch).")
+        comp2off = {"X": 0, "Y": 1, "Z": 2}
+
+        aktive: List[int] = []
+        for pt, comp in auswahl:
+            spt = str(pt)
+            c = comp.upper()
+            if spt not in punkt2pos:
+                raise KeyError(f"Punkt '{pt}' nicht in liste_punktnummern.")
+            if c not in comp2off:
+                raise ValueError(f"Komponente '{comp}' ungültig. Nur X,Y,Z.")
+            p = punkt2pos[spt]
+            aktive.append(3 * p + comp2off[c])
+
+        # Duplikate entfernen
+        out, seen = [], set()
+        for i in aktive:
+            if i not in seen:
+                seen.add(i)
+                out.append(i)
+        return out
+
+    @staticmethod
+    def auswahlmatrix_E(u: int, aktive_unbekannte_indices: Iterable[int]) -> sp.Matrix:
+        E = sp.zeros(u, u)
+        for idx in aktive_unbekannte_indices:
+            i = int(idx)
+            if not (0 <= i < u):
+                raise IndexError(f"Aktiver Index {i} außerhalb [0,{u-1}]")
+            E[i, i] = 1
+        return E
+
+    @staticmethod
+    def raenderungsmatrix_G(
+        x0: sp.Matrix,
+        liste_punktnummern: Sequence[str],
+        *,
+        mit_massstab: bool = True,
+        layout: str = "XYZ",
    ) -> sp.Matrix:
+        if x0.cols != 1:
+            raise ValueError("x0 muss Spaltenvektor sein.")
+        layout = layout.upper()
+        if layout != "XYZ":
+            raise ValueError("Nur layout='XYZ' unterstützt (wie bei euch).")

-    u = x0.rows
-    d = 7 if mit_massstab else 6
-    G = sp.zeros(u, d)
+        nP = len(liste_punktnummern)
+        u = x0.rows
+        d = 7 if mit_massstab else 6
+        G = sp.zeros(u, d)

-    # --- Translationen ---
-    for i in idx_X:
-        G[i, 0] = 1
-    for i in idx_Y:
-        G[i, 1] = 1
-    for i in idx_Z:
-        G[i, 2] = 1
+        for p in range(nP):
+            ix, iy, iz = 3*p, 3*p+1, 3*p+2
+            xi, yi, zi = x0[ix, 0], x0[iy, 0], x0[iz, 0]

-    # --- Rotationen ---
-    # Rotation um X-Achse
-    for iy, iz in zip(idx_Y, idx_Z):
-        zi = x0[iz, 0]
-        yi = x0[iy, 0]
-        G[iy, 3] = -zi
-        G[iz, 3] = yi
+            # Translationen
+            G[ix, 0] = 1
+            G[iy, 1] = 1
+            G[iz, 2] = 1

-    # Rotation um Y-Achse
-    for ix, iz in zip(idx_X, idx_Z):
-        zi = x0[iz, 0]
-        xi = x0[ix, 0]
-        G[ix, 4] = zi
-        G[iz, 4] = -xi
+            # Rotationen
+            G[iy, 3] = -zi; G[iz, 3] =  yi     # Rx
+            G[ix, 4] =  zi; G[iz, 4] = -xi     # Ry
+            G[ix, 5] = -yi; G[iy, 5] =  xi     # Rz

-    # Rotation um Z-Achse
-    for ix, iy in zip(idx_X, idx_Y):
-        yi = x0[iy, 0]
-        xi = x0[ix, 0]
-        G[ix, 5] = -yi
-        G[iy, 5] = xi
+            # Maßstab
+            if mit_massstab:
+                G[ix, 6] = xi
+                G[iy, 6] = yi
+                G[iz, 6] = zi

-    # --- Maßstab ---
-    if mit_massstab:
-        for ix, iy, iz in zip(idx_X, idx_Y, idx_Z):
-            xi = x0[ix, 0]
-            yi = x0[iy, 0]
-            zi = x0[iz, 0]
-            G[ix, 6] = xi
-            G[iy, 6] = yi
-            G[iz, 6] = zi
-    return G
+        return G

+    @staticmethod
+    def berechne_dx_geraendert(N: sp.Matrix, n: sp.Matrix, Gi: sp.Matrix) -> sp.Matrix:
+        if N.rows != N.cols:
+            raise ValueError("N muss quadratisch sein.")
+        if n.cols != 1:
+            raise ValueError("n muss Spaltenvektor sein.")
+        if Gi.rows != N.rows:
+            raise ValueError("Gi hat falsche Zeilenzahl.")

-def auswahlmatrix_E(u: int, aktive_unbekannte_indices: Iterable[int]) -> sp.Matrix:
-    E = sp.zeros(u, u)
-    for idx in aktive_unbekannte_indices:
-        E[int(idx), int(idx)] = 1
-    return E
+        u = N.rows
+        d = Gi.cols
+        K = N.row_join(Gi)
+        K = K.col_join(Gi.T.row_join(sp.zeros(d, d)))
+        rhs = n.col_join(sp.zeros(d, 1))
+        sol = K.LUsolve(rhs)
+        return sol[:u, :]

+    @staticmethod
+    def weiches_datum(
+        A: sp.Matrix,
+        dl: sp.Matrix,
+        Q_ll: sp.Matrix,
+        x0: sp.Matrix,
+        anschluss_indices: Sequence[int],
+        anschluss_werte: sp.Matrix,
+        Sigma_AA: Optional[sp.Matrix] = None,
+    ) -> Tuple[sp.Matrix, sp.Matrix, sp.Matrix]:
+        if dl.cols != 1 or x0.cols != 1:
+            raise ValueError("dl und x0 müssen Spaltenvektoren sein.")
+        if A.rows != dl.rows:
+            raise ValueError("A.rows muss dl.rows entsprechen.")
+        if A.cols != x0.rows:
+            raise ValueError("A.cols muss x0.rows entsprechen.")
+        if Q_ll.rows != Q_ll.cols or Q_ll.rows != A.rows:
+            raise ValueError("Q_ll muss (n×n) sein und zu A.rows passen.")

-def teilspurminimierung_Gi(G: sp.Matrix, E: sp.Matrix) -> sp.Matrix:
-    Gi = E * G
-    return Gi
+        u = A.cols
+        idx = [int(i) for i in anschluss_indices]
+        m = len(idx)

+        if anschluss_werte.cols != 1 or anschluss_werte.rows != m:
+            raise ValueError("anschluss_werte muss (m×1) sein.")
+        if Sigma_AA is None:
+            Sigma_AA = sp.eye(m)
+        if Sigma_AA.rows != m or Sigma_AA.cols != m:
+            raise ValueError("Sigma_AA muss (m×m) sein.")

-def berechne_dx_geraendert(N: sp.Matrix, n: sp.Matrix, Gi: sp.Matrix) -> sp.Matrix:
-    u = N.rows
-    d = Gi.shape[1]
+        A_A = sp.zeros(m, u)
+        for r, j in enumerate(idx):
+            if not (0 <= j < u):
+                raise IndexError(f"Anschluss-Index {j} außerhalb [0,{u-1}]")
+            A_A[r, j] = 1

-    K = N.row_join(Gi)
-    K = K.col_join(Gi.T.row_join(sp.zeros(d, d)))
+        x0_A = sp.Matrix([[x0[j, 0]] for j in idx])
+        dl_A = anschluss_werte - x0_A

-    rhs = n.col_join(sp.zeros(d, 1))
+        A_ext = A.col_join(A_A)
+        dl_ext = dl.col_join(dl_A)

-    sol = K.LUsolve(rhs)
-    dx = sol[:u, :]
-    return dx
+        Q_ext = sp.zeros(Q_ll.rows + m, Q_ll.cols + m)
+        Q_ext[:Q_ll.rows, :Q_ll.cols] = Q_ll
+        Q_ext[Q_ll.rows:, Q_ll.cols:] = Sigma_AA
+
+        return A_ext, dl_ext, Q_ext
--- a/Parameterschaetzung.py
+++ b/Parameterschaetzung.py
@@ -1,29 +1,62 @@
-from Datumsfestlegung import *
 from Stochastisches_Modell import StochastischesModell
+from Netzqualität_Genauigkeit import Genauigkeitsmaße
+from Datumsfestlegung import Datumsfestlegung
 import sympy as sp
 import Export
-import Netzqualität_Genauigkeit


-def ausgleichung_global(A, dl, stoch_modell: StochastischesModell):
+def ausgleichung_global(
+    A: sp.Matrix,
+    dl: sp.Matrix,
+    Q_ll: sp.Matrix,
+    x0: sp.Matrix,
+    idx_X, idx_Y, idx_Z,
+    anschluss_indices,
+    anschluss_werte,
+    Sigma_AA,
+):
+    # 1) Datumsfestlegung (weiches Datum) System erweitern
+    A_ext, dl_ext, Q_ext = Datumsfestlegung.weiches_datum(
+        A=A,
+        dl=dl,
+        Q_ll=Q_ll,
+        x0=x0,
+        anschluss_indices=anschluss_indices,
+        anschluss_werte=anschluss_werte,
+        Sigma_AA=Sigma_AA,
+    )

-    #Q_ll, P = stoch_modell.berechne_Qll_P()        #Kofaktormatrix und P-Matrix
-    P = sp.eye(A.shape[0])
-    N = A.T * P * A                                 #Normalgleichungsmatrix N
-    Q_xx = N.inv()                                  #Kofaktormatrix der Unbekannten Qxx
-    n = A.T * P * dl                                #Absolutgliedvektor n
+    # 2) Gewichtsmatrix P
+    P = StochastischesModell.berechne_P(Q_ext)

-    dx = N.LUsolve(n)                               #Zuschlagsvektor dx
+    # 3) Normalgleichungsmatrix N und Absolutgliedvektor n
+    N = A_ext.T * P * A_ext
+    n = A_ext.T * P * dl_ext

-    v = dl - A * dx                                 #Residuenvektor v
+    # 4) Zuschlagsvektor dx
+    dx = N.LUsolve(n)

+    # 5) Residuenvektor v
+    v = dl - A * dx
+
+    # 6) Kofaktormatrix der Unbekannten Q_xx
+    Q_xx = StochastischesModell.berechne_Q_xx(N)
+
+    # 7) Kofaktormatrix der Beobachtungen Q_ll_dach
    Q_ll_dach = A * Q_xx * A.T
-    Q_vv = stoch_modell.berechne_Qvv(A, P, Q_xx)          #Kofaktormatrix der Verbesserungen Qvv
-    R = stoch_modell.berechne_R(Q_vv, P)                  #Redundanzmatrix R
-    r = stoch_modell.berechne_r(R)                        #Redundanzanteile als Vektor r
-    redundanzanteile = A.shape[0] - A.shape[1] #n-u+d
-    soaposteriori = Netzqualität_Genauigkeit.Genauigkeitsmaße.s0apost(v, P, redundanzanteile)

+    # 8) Kofaktormatrix der Verbesserungen Q_vv
+    Q_vv = StochastischesModell.berechne_Qvv(A, P, Q_xx)
+
+    # 9) Redundanzmatrix R und Redundanzanteile r
+    R = StochastischesModell.berechne_R(Q_vv, P)                  #Redundanzmatrix R
+    r = StochastischesModell.berechne_r(R)                        #Redundanzanteile als Vektor r
+    redundanzanteile = A.shape[0] - A.shape[1] #n-u+d
+
+    # 10) s0 a posteriori
+    soaposteriori = Genauigkeitsmaße.s0apost(v, P, redundanzanteile)
+
+    # 11) Ausgabe
    dict_ausgleichung = {
        "dx": dx,
        "v": v,
@@ -38,68 +71,65 @@ def ausgleichung_global(A, dl, stoch_modell: StochastischesModell):
    }

    Export.Export.ausgleichung_to_datei(r"Zwischenergebnisse\Ausgleichung_Iteration0.csv", dict_ausgleichung)
-
    return dict_ausgleichung, dx



+
 def ausgleichung_lokal(
    A: sp.Matrix,
    dl: sp.Matrix,
-    stoch_modell: StochastischesModell,
+    Q_ll: sp.Matrix,
    x0: sp.Matrix,
    idx_X, idx_Y, idx_Z,
    aktive_unbekannte_indices,
    mit_massstab: bool = True,
 ):
-    # 1) Gewichte
-    Q_ll, P = stoch_modell.berechne_Qll_P()
-    # Debug-Option:
-    # P = sp.eye(A.rows)
+    # 1) Gewichtsmatrix P
+    P = StochastischesModell.berechne_P(Q_ll)

-    # 2) Normalgleichungen
+    # 2) Normalgleichungsmatrix N und Absolutgliedvektor n
    N = A.T * P * A
    n = A.T * P * dl

-    # 3) Datum (G, E, Gi)
-    G = raenderungsmatrix_G(x0, idx_X, idx_Y, idx_Z, mit_massstab=mit_massstab)
-    E = auswahlmatrix_E(u=A.cols, aktive_unbekannte_indices=aktive_unbekannte_indices)
+    # 3) Datumsfestlegung (Teilspurminimierung)
+    G = Datumsfestlegung.raenderungsmatrix_G(x0, liste_punktnummern, mit_massstab=mit_massstab)
+    aktive = Datumsfestlegung.datumskomponenten(auswahl, liste_punktnummern)
+    E = Datumsfestlegung.auswahlmatrix_E(u=A.cols, aktive_unbekannte_indices=aktive)
    Gi = E * G

-    # 4) Geränderte Lösung (dx)
-    dx = berechne_dx_geraendert(N, n, Gi)
+    # 4) Zuschlagsvektor dx
+    dx = Datumsfestlegung.berechne_dx_geraendert(N, n, Gi)

-    # 5) Residuen
+    # 5) Residuenvektor v
    v = dl - A * dx

-    # 6) KORREKTE Q_xx für gerändertes Problem:
-    # Q_xx = N^{-1} - N^{-1}Gi (Gi^T N^{-1} Gi)^{-1} Gi^T N^{-1}
-    # numerisch besser via LUsolve statt inv:
-    N_inv = N.inv()  # wenn N groß ist, kann man das unten auch ohne inv machen (siehe Hinweis)
+    # 6) Kofaktormatrix der Unbekannten Q_xx
+    N_inv = N.inv()
    N_inv_G = N_inv * Gi
    S = Gi.T * N_inv_G
    S_inv = S.inv()
    Q_xx = N_inv - N_inv_G * S_inv * N_inv_G.T

-    # 7) Q_lhat_lhat und Q_vv
+    # 7) Kofaktormatrix der Beobachtungen Q_ll_dach
    Q_lhat_lhat = A * Q_xx * A.T
+
+    # 8) Kofaktormatrix der Verbesserungen Q_vv
    Q_vv = P.inv() - Q_lhat_lhat

-    # 8) Redundanzmatrix und -anteile
+    # 9) Redundanzmatrix R, Redundanzanteile r, Redundanz
    R = Q_vv * P
    r_vec = sp.Matrix(R.diagonal())
-
-    # 9) Freiheitsgrade (Redundanz gesamt)
    n_beob = A.rows
    u = A.cols
    d = Gi.shape[1]
    r_gesamt = n_beob - u + d

-    # 10) sigma0 a posteriori
-    omega = float((v.T * P * v)[0, 0])
-    sigma0_hat = (omega / float(r_gesamt)) ** 0.5
+    # 10) s0 a posteriori
+    sigma0_apost = Genauigkeitsmaße.s0apost(v, P, r_gesamt)

-    return {
+    # 11) Ausgabe
+    dict_ausgleichung_lokal = {
        "dx": dx,
        "v": v,
        "Q_ll": Q_ll,
@@ -111,7 +141,10 @@ def ausgleichung_lokal(
        "R": R,
        "r": r_vec,
        "r_gesamt": r_gesamt,
-        "sigma0_hat": sigma0_hat,
+        "sigma0_apost": sigma0_apost,
        "G": G,
        "Gi": Gi,
-    }
+    }
+
+    Export.Export.ausgleichung_to_datei(r"Zwischenergebnisse\Ausgleichung_Iteration0_lokal.csv", dict_ausgleichung_lokal)
+    return dict_ausgleichung_lokal, dx
--- a/Stochastisches_Modell.py
+++ b/Stochastisches_Modell.py
@@ -49,11 +49,18 @@ class StochastischesModell:
        return Q_ll


-    def berechne_P(Q_ll):
+    def berechne_P(Q_ll: sp.Matrix) -> sp.Matrix:
        P = Q_ll.inv()
        return P


+    def berechne_Q_xx(N: sp.Matrix) -> sp.Matrix:
+        if N.rows != N.cols:
+            raise ValueError("N muss eine quadratische Matrix sein")
+        Q_xx = N.inv()
+        return Q_xx
+
+
    def berechne_Qvv(self, A: sp.Matrix, P: sp.Matrix, Q_xx: sp.Matrix) -> sp.Matrix:
        Q_vv = P.inv() - A * Q_xx * A.T
        return Q_vv                                 #Kofaktormatrix der Beobachtungsresiduen