diff --git a/src/site/resources/userguide/TrajectoryDeterminationProblem.java b/src/site/resources/userguide/TrajectoryDeterminationProblem.java
new file mode 100644
index 000000000..32a0975f1
--- /dev/null
+++ b/src/site/resources/userguide/TrajectoryDeterminationProblem.java
@@ -0,0 +1,187 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+import org.apache.commons.math.estimation.EstimationException;
+import org.apache.commons.math.estimation.EstimatedParameter;
+import org.apache.commons.math.estimation.EstimationProblem;
+import org.apache.commons.math.estimation.LevenbergMarquardtEstimator;
+import org.apache.commons.math.estimation.SimpleEstimationProblem;
+import org.apache.commons.math.estimation.WeightedMeasurement;
+
+public class TrajectoryDeterminationProblem extends SimpleEstimationProblem {
+
+ public static void main(String[] args) {
+ try {
+ TrajectoryDeterminationProblem problem =
+ new TrajectoryDeterminationProblem(0.0, 100.0, 800.0, 1.0, 0.0);
+
+ double[][] distances = {
+ { 0.0, 806.5849 }, { 20.0, 796.8148 }, { 40.0, 791.0833 }, { 60.0, 789.6712 },
+ { 80.0, 793.1334 }, { 100.0, 797.7248 }, { 120.0, 803.2785 }, { 140.0, 813.4939 },
+ { 160.0, 826.9295 }, { 180.0, 844.0640 }, { 200.0, 863.3829 }, { 220.0, 883.3143 },
+ { 240.0, 908.6867 }, { 260.0, 934.8561 }, { 280.0, 964.0730 }, { 300.0, 992.1033 },
+ { 320.0, 1023.998 }, { 340.0, 1057.439 }, { 360.0, 1091.912 }, { 380.0, 1125.968 },
+ { 400.0, 1162.789 }, { 420.0, 1201.517 }, { 440.0, 1239.176 }, { 460.0, 1279.347 } };
+ for (int i = 0; i < distances.length; ++i) {
+ problem.addDistanceMeasurement(1.0, distances[i][0], distances[i][1]);
+ };
+
+ double[][] angles = {
+ { 10.0, 1.415423 }, { 30.0, 1.352643 }, { 50.0, 1.289290 }, { 70.0, 1.225249 },
+ { 90.0, 1.161203 }, {110.0, 1.098538 }, {130.0, 1.036263 }, {150.0, 0.976052 },
+ {170.0, 0.917921 }, {190.0, 0.861830 }, {210.0, 0.808237 }, {230.0, 0.757043 },
+ {250.0, 0.708650 }, {270.0, 0.662949 }, {290.0, 0.619903 }, {310.0, 0.579160 },
+ {330.0, 0.541033 }, {350.0, 0.505590 }, {370.0, 0.471746 }, {390.0, 0.440155 },
+ {410.0, 0.410522 }, {430.0, 0.382701 }, {450.0, 0.356957 }, {470.0, 0.332400 } };
+ for (int i = 0; i < angles.length; ++i) {
+ problem.addAngularMeasurement(3.0e7, angles[i][0], angles[i][1]);
+ };
+
+ LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
+ estimator.estimate(problem);
+ System.out.println("initial position: " + problem.getX0() + " " + problem.getY0());
+ System.out.println("velocity: " + problem.getVx0() + " " + problem.getVy0());
+
+ } catch (EstimationException ee) {
+ System.err.println(ee.getMessage());
+ }
+ }
+
+ public TrajectoryDeterminationProblem(double t0,
+ double x0Guess, double y0Guess,
+ double vx0Guess, double vy0Guess) {
+ this.t0 = t0;
+ x0 = new EstimatedParameter( "x0", x0Guess);
+ y0 = new EstimatedParameter( "y0", y0Guess);
+ vx0 = new EstimatedParameter("vx0", vx0Guess);
+ vy0 = new EstimatedParameter("vy0", vy0Guess);
+
+ // inform the base class about the parameters
+ addParameter(x0);
+ addParameter(y0);
+ addParameter(vx0);
+ addParameter(vy0);
+
+ }
+
+ public double getX0() {
+ return x0.getEstimate();
+ }
+
+ public double getY0() {
+ return y0.getEstimate();
+ }
+
+ public double getVx0() {
+ return vx0.getEstimate();
+ }
+
+ public double getVy0() {
+ return vy0.getEstimate();
+ }
+
+ public void addAngularMeasurement(double wi, double ti, double ai) {
+ // let the base class handle the measurement
+ addMeasurement(new AngularMeasurement(wi, ti, ai));
+ }
+
+ public void addDistanceMeasurement(double wi, double ti, double di) {
+ // let the base class handle the measurement
+ addMeasurement(new DistanceMeasurement(wi, ti, di));
+ }
+
+ public double x(double t) {
+ return x0.getEstimate() + (t - t0) * vx0.getEstimate();
+ }
+
+ public double y(double t) {
+ return y0.getEstimate() + (t - t0) * vy0.getEstimate();
+ }
+
+ private class AngularMeasurement extends WeightedMeasurement {
+
+ public AngularMeasurement(double weight, double t, double angle) {
+ super(weight, angle);
+ this.t = t;
+ }
+
+ public double getTheoreticalValue() {
+ return Math.atan2(y(t), x(t));
+ }
+
+ public double getPartial(EstimatedParameter parameter) {
+ double xt = x(t);
+ double yt = y(t);
+ double r = Math.sqrt(xt * xt + yt * yt);
+ double u = yt / (r + xt);
+ double c = 2 * u / (1 + u * u);
+ if (parameter == x0) {
+ return -c;
+ } else if (parameter == vx0) {
+ return -c * t;
+ } else if (parameter == y0) {
+ return c * xt / yt;
+ } else {
+ return c * t * xt / yt;
+ }
+ }
+
+ private final double t;
+ private static final long serialVersionUID = -5990040582592763282L;
+
+ }
+
+ private class DistanceMeasurement extends WeightedMeasurement {
+
+ public DistanceMeasurement(double weight, double t, double angle) {
+ super(weight, angle);
+ this.t = t;
+ }
+
+ public double getTheoreticalValue() {
+ double xt = x(t);
+ double yt = y(t);
+ return Math.sqrt(xt * xt + yt * yt);
+ }
+
+ public double getPartial(EstimatedParameter parameter) {
+ double xt = x(t);
+ double yt = y(t);
+ double r = Math.sqrt(xt * xt + yt * yt);
+ if (parameter == x0) {
+ return xt / r;
+ } else if (parameter == vx0) {
+ return xt * t / r;
+ } else if (parameter == y0) {
+ return yt / r;
+ } else {
+ return yt * t / r;
+ }
+ }
+
+ private final double t;
+ private static final long serialVersionUID = 3257286197740459503L;
+
+ }
+
+ private double t0;
+ private EstimatedParameter x0;
+ private EstimatedParameter y0;
+ private EstimatedParameter vx0;
+ private EstimatedParameter vy0;
+
+}
diff --git a/src/site/resources/userguide/estimation-class-diagram.png b/src/site/resources/userguide/estimation-class-diagram.png
new file mode 100644
index 000000000..fe9b6786a
Binary files /dev/null and b/src/site/resources/userguide/estimation-class-diagram.png differ
diff --git a/src/site/resources/userguide/estimation-sequence-diagram.png b/src/site/resources/userguide/estimation-sequence-diagram.png
new file mode 100644
index 000000000..f678bb283
Binary files /dev/null and b/src/site/resources/userguide/estimation-sequence-diagram.png differ
diff --git a/xdocs/userguide/estimation.xml b/xdocs/userguide/estimation.xml
index c214ab019..9dadf7368 100644
--- a/xdocs/userguide/estimation.xml
+++ b/xdocs/userguide/estimation.xml
@@ -1,22 +1,22 @@
+ http://www.apache.org/licenses/LICENSE-2.0
+
+ Unless required by applicable law or agreed to in writing, software
+ distributed under the License is distributed on an "AS IS" BASIS,
+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ See the License for the specific language governing permissions and
+ limitations under the License.
+-->
+
- The estimation package provides classes to fit some non-linear model
- to available observations depending on it. These problems are commonly
- called estimation problems.
+ The estimation package provides classes to fit some non-linear
+ model to available observations depending on it. These
+ problems are commonly called estimation problems.
- The estimation problems considered here are parametric problems where
- a user-provided model depends on initially unknown scalar parameters and
- several measurements made on values that depend on the model are available.
- As examples, one can consider the center and radius of a circle given
- points approximately lying on a ring, or a satellite orbit given range,
- range-rate and angular measurements from various ground stations.
+ The estimation problems considered here are parametric
+ problems where a user-provided model depends on initially
+ unknown scalar parameters and several measurements made on
+ values that depend on the model are available. As examples,
+ one can consider the center and radius of a circle given
+ points approximately lying on a ring, or a satellite orbit
+ given range, range-rate and angular measurements from various
+ ground stations.
- One important class of estimation problems is weighted least squares problems.
- They basically consist in finding the values for some parameters pk
- such that a cost function J =
-
- sum(wi ri2)
- is minimized. The various ri terms represent the deviation
- ri = mesi - modi between the measurements and
- the parameterized models. The wi factors are the measurements weights,
- they are often chosen either all equal to 1.0 or proportional to the inverse of
- the variance of the measurement type. The solver adjusts the values of the
- estimated parameters pk which are not bound (i.e. the free parameters).
- It does not touch the parameters which have been put in a bound state by the user.
+ One important class of estimation problems is weighted least
+ squares problems. They basically consist in finding the values
+ for some parameters pk such that a cost function
+ J = sum(wiri2) is minimized.
+ The various ri terms represent the deviation
+ ri = mesi - modi
+ between the measurements and the parameterized models. The
+ wi factors are the measurements weights, they are often
+ chosen either all equal to 1.0 or proportional to the inverse of the
+ variance of the measurement type. The solver adjusts the values of
+ the estimated parameters pk which are not bound (i.e. the
+ free parameters). It does not touch the parameters which have been
+ put in a bound state by the user.
- The aim of this package is similar to the aim of the optimization package, but the
- algorithms are entirely differents as:
+ The aim of this package is similar to the aim of the
+ optimization package, but the algorithms are entirely
+ different as:
- The - org.apache.commons.math.estimation.EstimatedParameter class to represent each - estimated parameter. The parameters are set up with a guessed value which will be - adjusted by the solver until a best fit is achieved. It is possible to change which - parameters are modified and which are preserved thanks to a bound property. Such - settings are often needed by expert users to handle contingency cases with very - low observability. + The problem modeling is the most important part for the + user. Understanding it is the key to proper use of the + package. One interface and two classes are provided for this + purpose: + EstimationProblem, + EstimatedParameter and + WeightedMeasurement.
-- The user extends the - org.apache.commons.math.estimation.WeightedMeasurement abstract class to define its - own measurements. Each measurement types should have its own implementing class, for - example in the satellite example above , the user should define three classes, one - for range measurements, one for range-rates measurements and one for angular measurements. - Each measurement would correspond to an instance of the appropriate class, set up with - the date, a reference to the ground station, the weight and the measured value. -
-- The package provides two common - org.apache.commons.math.estimation.Estimator implementations to solve weighted - least squares problems. The first one is based on the - Gauss-Newton method. - The second one is based on the - Levenberg-Marquardt - method. The first one is best suited when a good approximation of the parameters is known while the second one - is more robust and can handle starting points far from the solution. + Consider the following example problem: we want to determine the + linear trajectory of a sailing ship by performing angular and + distance measurements from an observing spot on the shore. The + problem model is represented by two equations:
-
+ x(t) = x0+(t-t0)vx0
+ y(t) = y0+(t-t0)vy0
+
+ These two equations depend on four parameters (x0, y0, + vx0 and vy0). We want to determine these four parameters. +
++ Assuming the observing spot is located at the origin of the coordinates + system and that the angular measurements correspond to the angle between + the x axis and the line of sight, the theoretical values of the angular + measurements at ti and of the distance measurements at + tj are modeled as follows: +
+
+ anglei,theo = atan2(y(ti), x(ti))
+ distancej,theo = sqrt(x(tj)2+y(tj)2)
+
+ The real observations generate a set of measurements values anglei,meas + and distancej,meas. +
++ The following class diagram shows one way to solve this problem using the + estimation package. The grey elements are already provided by the package + whereas the purple elements are developed by the user. +
+ +
+ The TrajectoryDeterminationProblem
class holds the linear model
+ equations x(t) and y(t). It delegate storage of the four parameters x0,
+ y0, vx0 and vy0 and of the various measurements
+ anglei,meas and distancej,meas to its base class
+ SimpleEstimationProblem
. Since the theoretical values of the measurements
+ anglei,theo and distancej,theo depend on the linear model,
+ the two classes AngularMeasurement
and DistanceMeasurement
+ are implemented as internal classes, thus having access to the equations of the
+ linear model and to the parameters.
+
+ Here are the various parts of the TrajectoryDeterminationProblem.java
+ source file. This example, with an additional main
method is
+ available here.
+
serialVersionUID
static fields are present because
+ the WeightedMeasurement
class implements the
+ Serializable
interface.
+
+
+
+ Solving the problem is simply a matter of choosing an implementation
+ of the
+ Estimator interface and to pass the problem instance to its estimate
+ method. Two implementations are already provided by the library:
+ GaussNewtonEstimator and
+ LevenbergMarquardtEstimator. The first one implements a simple Gauss-Newton
+ algorithm, which is sufficient when the starting point (initial guess) is close
+ enough to the solution. The second one implements a more complex Levenberg-Marquardt
+ algorithm which is more robust when the initial guess is far from the solution.
+
+ The following sequence diagram explains roughly what occurs under the hood
+ in the estimate
method.
+
+ Basically, the estimator first retrieves the parameters and the measurements. + The estimation loop is based on the gradient of the sum of the squares of the + residuals, hence, the estimators get the various partial derivatives of all + measurements with respect to all parameters. A new state hopefully globally + reducing the residuals is estimated, and the parameters value are updated. + This estimation loops stops when either the convergence conditions are met + or the maximal number of iterations is exceeded. +
++ One important tuning parameter for weighted least-squares solving is the + weight attributed to each measurement. This weights has two purposes: +
++ The weight is a multiplicative factor for the square of the residuals. + A common choice is to use the inverse of the variance of the measurements error + as the weighting factor for all measurements for one type. On our sailing ship + example, we may have a range measurements accuracy of about 1 meter and an angular + measurements accuracy of about 0.01 degree, or 1.7 10-4 radians. So we + would use w=1.0 for distance measurements weight and w=3 107 for + angular measurements weight. If we knew that the measurements quality is bad + at tracking start because of measurement system warm-up delay for example, then + we would reduce the weight for the first measurements and use for example + w=0.1 and w=3 106 respectively, depending on the type. +
+
+ After a problem has been set up, it is possible to fine tune the
+ way it will be solved. For example, it may appear the measurements are not
+ sufficient to get some parameters with sufficient confidence due to observability
+ problems. It is possible to fix some parameters in order to prevent the solver
+ from changing them. This is realized by passing true
to the
+ setBound
method of the parameter.
+
+ It is also possible to ignore some measurements by passing true
to the
+ setIgnored
method of the measurement. A typical use is to
+