diff --git a/.idea/Masterprojekt-Campusnetz.iml b/.idea/Masterprojekt-Campusnetz.iml index 1d2fcdf..89b2bd1 100644 --- a/.idea/Masterprojekt-Campusnetz.iml +++ b/.idea/Masterprojekt-Campusnetz.iml @@ -4,7 +4,7 @@ - + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml index ba45cb8..590a59e 100644 --- a/.idea/misc.xml +++ b/.idea/misc.xml @@ -3,5 +3,5 @@ - + \ No newline at end of file diff --git a/Parameterschaetzung.py b/Parameterschaetzung.py index 28f4d1b..e555319 100644 --- a/Parameterschaetzung.py +++ b/Parameterschaetzung.py @@ -1,77 +1,29 @@ -from typing import Dict, Any -import sympy as sp from Stochastisches_Modell import StochastischesModell -def iterative_ausgleichung( - A: sp.Matrix, - l: sp.Matrix, - modell: StochastischesModell, - max_iter: int = 100, - tol: float = 1e-3, -) -> Dict[str, Any]: +def ausgleichung(A, dl, stoch_modell: StochastischesModell, P): - ergebnisse_iter = [] #Liste für Zwischenergebnisse + Q_ll, P = stoch_modell.berechne_Qll_P() #Kofaktormatrix und P-Matrix + N = A.T * P * A #Normalgleichungsmatrix N + Q_xx = N.inv() #Kofaktormatrix der Unbekannten Qxx + n = A.T * P * dl #Absolutgliedvektor n - for it in range(max_iter): - Q_ll, P = modell.berechne_Qll_P() #Stochastisches Modell: Qll und P berechnen + dx = N.LUsolve(n) #Zuschlagsvektor dx - N = A.T * P * A #Normalgleichungsmatrix N - Q_xx = N.inv() #Kofaktormatrix der Unbekannten Qxx - n = A.T * P * l #Absolutgliedvektor n + v = dl - A * dx #Residuenvektor v - dx = N.LUsolve(n) #Zuschlagsvektor dx - - v = l - A * dx #Residuenvektor v - - Q_vv = modell.berechne_Qvv(A, P, Q_xx) #Kofaktormatrix der Verbesserungen Qvv - R = modell.berechne_R(Q_vv, P) #Redundanzmatrix R - r = modell.berechne_r(R) #Redundanzanteile als Vektor r - - sigma_hat = modell.berechne_vks(v, P, r) #Varianzkomponentenschätzung durchführen - - ergebnisse_iter.append({ #Zwischenergebnisse speichern in Liste - "iter": it + 1, - "Q_ll": Q_ll, - "P": P, - "N": N, - "Q_xx": Q_xx, - "dx": dx, - "v": v, - "Q_vv": Q_vv, - "R": R, - "r": r, - "sigma_hat": sigma_hat, - "sigma0_groups": dict(modell.sigma0_groups), - }) - - # --- Abbruchkriterium --- - if sigma_hat: - max_rel_change = 0.0 - for g, new_val in sigma_hat.items(): - old_val = modell.sigma0_groups.get(g, 1.0) - if old_val != 0: - rel = abs(new_val - old_val) / abs(old_val) - max_rel_change = max(max_rel_change, rel) - - if max_rel_change < tol: - print(f"Konvergenz nach {it + 1} Iterationen erreicht (max. rel. Änderung = {max_rel_change:.2e}).") - modell.update_sigma0_von_vks(sigma_hat) - break - - # Varianzfaktoren für nächste Iteration übernehmen - modell.update_sigma0_von_vks(sigma_hat) + 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 return { "dx": dx, "v": v, - "Q_ll": Q_ll, "P": P, "N": N, "Q_xx": Q_xx, + "Q_ll_dach": Q_ll_dach, "Q_vv": Q_vv, "R": R, "r": r, - "sigma_hat": sigma_hat, - "sigma0_groups": dict(modell.sigma0_groups), - "history": ergebnisse_iter, } \ No newline at end of file diff --git a/Stochastisches_Modell.py b/Stochastisches_Modell.py index 8b127d8..fc332b5 100644 --- a/Stochastisches_Modell.py +++ b/Stochastisches_Modell.py @@ -4,22 +4,22 @@ from typing import Dict, Tuple, Iterable @dataclass class StochastischesModell: - sigma_beob: Iterable[float] #σ der einzelnen Beobachtung - group_beob: Iterable[int] #Gruppenzugehörigkeit jeder Beobachtung (Distanz, Richtung, GNSS, Nivellement,...,) - sigma0_groups: Dict[int, float] = field(default_factory=dict) #σ0² für jede Gruppe + sigma_beob: Iterable[float] #σ a priori der einzelnen Beobachtung + gruppe_beob: Iterable[int] #Gruppenzugehörigkeit jeder Beobachtung (Distanz, Richtung, GNSS, Nivellement,...,) + sigma0_gruppe: Dict[int, float] = field(default_factory=dict) #σ0² für jede Gruppe def __post_init__(self): - self.sigma_beob = sp.Matrix(list(self.sigma_beob)) #Spaltenvektor - self.group_beob = sp.Matrix(list(self.group_beob)) #Spaltenvektor + self.sigma_beob = sp.Matrix(list(self.sigma_beob)) #Spaltenvektor + self.gruppe_beob = sp.Matrix(list(self.gruppe_beob)) #Spaltenvektor - if self.sigma_beob.rows != self.group_beob.rows: + if self.sigma_beob.rows != self.gruppe_beob.rows: raise ValueError("sigma_obs und group_ids müssen gleich viele Einträge haben.") - unique_groups = sorted({int(g) for g in self.group_beob}) #jede Beobachtungsgruppe wird genau einmal berücksichtigt + unique_groups = sorted({int(g) for g in self.gruppe_beob}) #jede Beobachtungsgruppe wird genau einmal berücksichtigt for g in unique_groups: - if g not in self.sigma0_groups: #Fehlende Gruppen mit σ_0j^2 = 1.0 - self.sigma0_groups[g] = 1.0 + if g not in self.sigma0_gruppe: #Fehlende Gruppen mit σ_0j^2 = 1.0 + self.sigma0_gruppe[g] = 1.0 @property @@ -31,11 +31,10 @@ class StochastischesModell: n = self.n_beob Q_ll = sp.zeros(n, n) P = sp.zeros(n, n) - for i in range(self.n_beob): sigma_i = self.sigma_beob[i, 0] #σ-Wert der i-ten Beobachtung holen - g = int(self.group_beob[i, 0]) #Gruppenzugehörigkeit der Beobachtung bestimmen - sigma0_sq = self.sigma0_groups[g] #Den Varianzfaktor der Gruppe holen + g = int(self.gruppe_beob[i, 0]) #Gruppenzugehörigkeit der Beobachtung bestimmen + sigma0_sq = self.sigma0_gruppe[g] #Den Varianzfaktor der Gruppe holen q_ii = sigma_i**2 #σ² berechnen Q_ll[i, i] = q_ii #Diagonale P[i, i] = 1 / (sigma0_sq * q_ii) #durch VKS nicht mehr P=Qll^-1 @@ -57,30 +56,4 @@ class StochastischesModell: r = sp.zeros(n, 1) for i in range(n): r[i, 0] = R[i, i] - return r #Redundanzanteile - - - def berechne_vks(self,v: sp.Matrix, P: sp.Matrix, r: sp.Matrix) -> Dict[int, float]: - if v.rows != self.n_beob: - raise ValueError("v passt nicht zur Anzahl der Beobachtungen.") - gruppen = sorted({int(g) for g in self.group_beob}) - sigma_gruppen: Dict[int, float] = {} - for g in gruppen: - idx = [i for i in range(self.n_beob) - if int(self.group_beob[i, 0]) == g] - if not idx: - continue - - v_g = sp.Matrix([v[i, 0] for i in idx]) - P_g = sp.zeros(len(idx), len(idx)) - for k, i_beob in enumerate(idx): - P_g[k, k] = P[i_beob, i_beob] - r_g = sum(r[i_beob, 0] for i_beob in idx) - sigma_gruppe_g = (v_g.T * P_g * v_g)[0, 0] / r_g - sigma_gruppen[g] = float(sigma_gruppe_g) - return sigma_gruppen - - - def update_sigma0_von_vks(self, sigma_hat: Dict[int, float]) -> None: - for g, val in sigma_hat.items(): - self.sigma0_groups[int(g)] = float(val) \ No newline at end of file + return r #Redundanzanteile \ No newline at end of file diff --git a/Tests_Michelle/Parameterschaetzung_müll.py b/Tests_Michelle/Parameterschaetzung_müll.py new file mode 100644 index 0000000..c9e03ec --- /dev/null +++ b/Tests_Michelle/Parameterschaetzung_müll.py @@ -0,0 +1,56 @@ +from typing import Dict, Any +import sympy as sp +from Stochastisches_Modell import StochastischesModell + +def iterative_ausgleichung( + A: sp.Matrix, + dl: sp.Matrix, + modell: StochastischesModell, +) -> Dict[str, Any]: + + ergebnisse_iter = [] #Liste für Zwischenergebnisse + + for it in range(max_iter): + Q_ll, P = modell.berechne_Qll_P() #Stochastisches Modell: Qll und P berechnen + + N = A.T * P * A #Normalgleichungsmatrix N + Q_xx = N.inv() #Kofaktormatrix der Unbekannten Qxx + n = A.T * P * dl #Absolutgliedvektor n + + dx = N.LUsolve(n) #Zuschlagsvektor dx + + v = dl - A * dx #Residuenvektor v + + Q_vv = modell.berechne_Qvv(A, P, Q_xx) #Kofaktormatrix der Verbesserungen Qvv + R = modell.berechne_R(Q_vv, P) #Redundanzmatrix R + r = modell.berechne_r(R) #Redundanzanteile als Vektor r + + ergebnisse_iter.append({ #Zwischenergebnisse speichern in Liste + "iter": it + 1, + "Q_ll": Q_ll, + "P": P, + "N": N, + "Q_xx": Q_xx, + "dx": dx, + "v": v, + "Q_vv": Q_vv, + "R": R, + "r": r, + "sigma_hat": sigma_hat, + "sigma0_groups": dict(modell.sigma0_groups), + }) + + return { + "dx": dx, + "v": v, + "Q_ll": Q_ll, + "P": P, + "N": N, + "Q_xx": Q_xx, + "Q_vv": Q_vv, + "R": R, + "r": r, + "sigma_hat": sigma_hat, + "sigma0_groups": dict(modell.sigma0_groups), + "history": ergebnisse_iter, + } \ No newline at end of file