diff --git a/Campusnetz.ipynb b/Campusnetz.ipynb index 70d32b0..c12849d 100644 --- 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": { diff --git a/Datumsfestlegung.py b/Datumsfestlegung.py index 8e703be..b4813e7 100644 --- 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 \ No newline at end of file + 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 diff --git a/Parameterschaetzung.py b/Parameterschaetzung.py index e1d5165..06d3ce5 100644 --- 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, - } \ No newline at end of file + } + + Export.Export.ausgleichung_to_datei(r"Zwischenergebnisse\Ausgleichung_Iteration0_lokal.csv", dict_ausgleichung_lokal) + return dict_ausgleichung_lokal, dx \ No newline at end of file diff --git a/Stochastisches_Modell.py b/Stochastisches_Modell.py index e30cefb..d544571 100644 --- 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