Small modifications to code; output

estimation/estimation_dynamical_models.ipynb

@@ -647,7 +647,7 @@
     "ax1.plot(time2plt_normalized, (mex_prop[:, 2] - mex_sim_obs[:, 2]), label=r'$\\Delta z$')\n",
     "ax1.plot(time2plt_normalized, np.linalg.norm((mex_prop[:, 0:3] - mex_sim_obs[:, 0:3]), axis=1), label=r'$||\\Delta X||$')\n",
     "\n",
-    "ax1.set_title(\"Element-wise difference with propagated observations\")\n",
+    "ax1.set_title(\"Element-wise difference between true and estimated states\")\n",
     "ax1.set_xlabel(r'$Time$ [days]')\n",
     "ax1.set_ylabel(r'$\\Delta X$ [m]')\n",
     "ax1.legend()\n",

estimation/estimation_dynamical_models.py

@@ -122,9 +122,9 @@
 
 ## Set Up the Observations
 """
+Having set the underlying environment model of the simulated orbit, we can define the observational model. This entails the addition all required ground stations, the definition of the observation links and types, as well as the precise simulation settings.
 """
 
-Having set the underlying environment model of the simulated orbit, we can define the observational model. This entails the addition all required ground stations, the definition of the observation links and types, as well as the precise simulation settings.
 ### Add a ground station
 """
 Following its real-world counterpart, our simulated Mars Express spacecraft will also be tracked using ESA's New Norcia (NNO) ESTRACK ground station. Located in North-East Australia, it will be set up with an altitude of 252m, 31.0482°S, 116.191°E.
@@ -433,7 +433,7 @@
 ax1.plot(time2plt_normalized, (mex_prop[:, 2] - mex_sim_obs[:, 2]), label=r'$\Delta z$')
 ax1.plot(time2plt_normalized, np.linalg.norm((mex_prop[:, 0:3] - mex_sim_obs[:, 0:3]), axis=1), label=r'$||\Delta X||$')
 
-ax1.set_title("Element-wise difference with propagated observations")
+ax1.set_title("Element-wise difference between true and estimated states")
 ax1.set_xlabel(r'$Time$ [days]')
 ax1.set_ylabel(r'$\Delta X$ [m]')
 ax1.legend()

estimation/full_estimation_example.ipynb

@@ -33,9 +33,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "b71677c3",
-   "metadata": {},
+   "metadata": {
+    "is_executing": true
+   },
    "outputs": [],
    "source": [
     "# Load required standard modules\n",
@@ -538,6 +540,13 @@
     "# Save the true parameters to later analyse the error\n",
     "truth_parameters = parameters_to_estimate.parameter_vector\n",
     "\n",
+    "# Perturb the initial state estimate from the truth (10 m in position; 0.1 m/s in velocity)\n",
+    "perturbed_parameters = truth_parameters.copy( )\n",
+    "for i in range(3):\n",
+    "    perturbed_parameters[i] += 10.0\n",
+    "    perturbed_parameters[i+3] += 0.01\n",
+    "parameters_to_estimate.parameter_vector = perturbed_parameters\n",
+    "\n",
     "# Create input object for the estimation\n",
     "convergence_checker = estimation.estimation_convergence_checker(maximum_iterations=4)\n",
     "estimation_input = estimation.EstimationInput(\n",
@@ -658,7 +667,7 @@
     "\n",
     "plt.figure(figsize=(9, 5))\n",
     "plt.title(\"Observations as a function of time\")\n",
-    "plt.scatter(observation_times / 3600.0, observations_list * constants.SPEED_OF_LIGHT)\n",
+    "plt.scatter(observation_times / 3600.0, observations_list)\n",
     "\n",
     "plt.xlabel(\"Time [hr]\")\n",
     "plt.ylabel(\"Range rate [m/s]\")\n",
@@ -697,16 +706,18 @@
    "source": [
     "residual_history = estimation_output.residual_history\n",
     "\n",
-    "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6), sharex=True, sharey=True)\n",
+    "fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6))\n",
     "subplots_list = [ax1, ax2, ax3, ax4]\n",
     "\n",
     "for i in range(4):\n",
     "    subplots_list[i].scatter(observation_times, residual_history[:, i])\n",
+    "    subplots_list[i].set_ylabel(\"Observation Residual [m/s]\")\n",
+    "    subplots_list[i].set_title(\"Iteration \"+str(i+1))\n",
+    "\n",
     "\n",
     "ax3.set_xlabel(\"Time since J2000 [s]\")\n",
     "ax4.set_xlabel(\"Time since J2000 [s]\")\n",
-    "ax1.set_ylabel(\"Observation Residual [m/s]\")\n",
-    "ax3.set_ylabel(\"Observation Residual [m/s]\")\n",
+    "\n",
     "\n",
     "plt.tight_layout()\n",
     "plt.show()"
@@ -762,6 +773,9 @@
     "\n",
     "plt.figure(figsize=(9,5))\n",
     "plt.hist(final_residuals, 25)\n",
+    "plt.xlabel('Final iteration range-rate residual [m/s]')\n",
+    "plt.ylabel('Occurences [-]')\n",
+    "plt.title('Histogram of residuals on final iteration')\n",
     "\n",
     "plt.tight_layout()\n",
     "plt.show()"

estimation/full_estimation_example.py

@@ -365,6 +365,14 @@
 # Save the true parameters to later analyse the error
 truth_parameters = parameters_to_estimate.parameter_vector
 
+# Perturb the initial state estimate from the truth (10 m in position; 0.1 m/s in velocity)
+perturbed_parameters = truth_parameters.copy( )
+for i in range(3):
+    perturbed_parameters[i] += 10.0
+    perturbed_parameters[i+3] += 0.01
+parameters_to_estimate.parameter_vector = perturbed_parameters
+
+
 # Create input object for the estimation
 convergence_checker = estimation.estimation_convergence_checker(maximum_iterations=4)
 estimation_input = estimation.EstimationInput(
@@ -412,7 +420,7 @@
 
 plt.figure(figsize=(9, 5))
 plt.title("Observations as a function of time")
-plt.scatter(observation_times / 3600.0, observations_list * constants.SPEED_OF_LIGHT)
+plt.scatter(observation_times / 3600.0, observations_list )
 
 plt.xlabel("Time [hr]")
 plt.ylabel("Range rate [m/s]")
@@ -429,16 +437,16 @@
 
 residual_history = estimation_output.residual_history
 
-fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6), sharex=True, sharey=True)
+fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(9, 6), sharex=True)
 subplots_list = [ax1, ax2, ax3, ax4]
 
 for i in range(4):
     subplots_list[i].scatter(observation_times, residual_history[:, i])
+    subplots_list[i].set_ylabel("Observation Residual [m/s]")
+    subplots_list[i].set_title("Iteration "+str(i+1))
 
 ax3.set_xlabel("Time since J2000 [s]")
 ax4.set_xlabel("Time since J2000 [s]")
-ax1.set_ylabel("Observation Residual [m/s]")
-ax3.set_ylabel("Observation Residual [m/s]")
 
 plt.tight_layout()
 plt.show()
@@ -459,7 +467,9 @@
 
 plt.figure(figsize=(9,5))
 plt.hist(final_residuals, 25)
-
+plt.xlabel('Final iteration range-rate residual [m/s]')
+plt.ylabel('Occurences [-]')
+plt.title('Histogram of residuals on final iteration')
 plt.tight_layout()
 plt.show()