MixedNoiseMeasurementModel

The MixedNoiseMeasurementModel interface describes the relationship between a system state xk∈RN and a noisy measurement yk∈RM in the form of yk=hk(xk,vk)+rk, with

  • measurement equation hk(xk,vk),

  • measurement noise vk∈RV, and

  • additive measurement noise rk∈RM.

The following example code can be found in the toolbox's examples.

Usage

In order to write a specific mixed noise measurement model, you have to

  • create a new subclass of MixedNoiseMeasurementModel,

  • implement hk(xk,vk) by implementing the abstract measurementEquation() method, and

  • set the measurement noise vk with the inherited setNoise() method.

  • set the additive measurement noise rk with the inherited setAdditiveNoise() method.

The noise has to be an instance of a Distribution subclass, e.g., a Gaussian. See also the list of available Probability Distributions.

You can access the set noise through the MixedNoiseMeasurementModel's properties.

Example

We implement a nonlinear measurement model for a 2D system state xk=[pk,qk]⊤. The measurement model is given by yk=hk(xk,vk)+rk=[pkq2kv2kp2k+3qkv3k]+rk, where vk is zero-mean white Gaussian noise with variance σ2=0.1 and rk zero-mean white Gaussian noise with covariance matrix R=I2.

The measurementEquation() method has to be implemented such that it can process an arbitrary number of L passed state samples [x(1)k,…,x(L)k] and corresponding noise samples [v(1)k,…,v(L)k] as it is done in the following. That is, measurementEquation() has to compute for each pair of state sample x(i)k the noise sample v(i)k the corresponding measurement y(i)k=hk(x(i)k,v(i)k) and return all measurements as matrix [y(1)k,…,y(L)k].

classdef MixedNoiseMeasModelExample < MixedNoiseMeasurementModel
    methods
        function obj = MixedNoiseMeasModelExample()
            % Specify the Gaussian system noise.
            % Of course, you do not have to call setNoise() and setAdditiveNoise() in the constructor.
            % You can also set/change the noise after creating a MixedNoiseMeasModelExample object.
            obj.setNoise(Gaussian(0, 0.1));
            obj.setAdditiveNoise(Gaussian(zeros(2, 1), eye(2)));
        end
        
        function measurements = measurementEquation(obj, stateSamples, noiseSamples)
            p = stateSamples(1, :);
            q = stateSamples(2, :);
            v = noiseSamples;
            
            measurements = [p .* q.^2 .* v.^2
                            p.^2 + 3*q .* v.^3];
        end
    end
end

If your measurement model has time-varying components, implement them as class properties and use these inside the measurementEquation(). Before performing a measurement update, you can adapt the measurement model by simply changing the class properties. For example, the measurement noise is handled in this way: by calling the setNoise() method before performing a measurement update, you can implement time-varying measurement noise.

Implement First-Order and Second-Order Derivatives of hk(xk,vk)

If you intend to use an estimator that relies on the first-order and maybe the second-order derivatives of the measurement equation hk(xk,vk), e.g., the EKF or the EKF2, you should consider implementing its analytical derivatives as well.

By default, the derivatives are automatically approximated in the derivative() method inherited from MixedNoiseMeasurementModel using the implemented measurementEquation() method and finite differences. However, the approximations might not be accurate enough or computationally expensive compared to an analytic implementation. In such a case, you can simply overwrite the derivative() method, which has to return the first-order and second-order derivatives for a provided nominal system state and nominal measurement noise, e.g., the prior state mean and the noise mean.

Example

We implement the first-order and second-order derivatives for the above example.

classdef MixedNoiseMeasModelExample < MixedNoiseMeasurementModel
    methods
        function obj = MixedNoiseMeasModelExample()
            % Specify the Gaussian system noise.
            % Of course, you do not have to call setNoise() and setAdditiveNoise() in the constructor.
            % You can also set/change the noise after creating a MixedNoiseMeasModelExample object.
            obj.setNoise(Gaussian(0, 0.1));
            obj.setAdditiveNoise(Gaussian(zeros(2, 1), eye(2)));
        end
        
        function measurements = measurementEquation(obj, stateSamples, noiseSamples)
            p = stateSamples(1, :);
            q = stateSamples(2, :);
            v = noiseSamples;
            
            measurements = [p .* q.^2 .* v.^2
                            p.^2 + 3*q .* v.^3];
        end
        
        function [stateJacobian, noiseJacobian, ...
                  stateHessians, noiseHessians] = derivative(obj, nominalState, nominalNoise)
            p = nominalState(1);
            q = nominalState(2);
            v = nominalNoise;
            
            % Jacobians
            stateJacobian = [q^2*v^2 2*p*q*v^2
                               2*p     3*v^3  ];
            
            noiseJacobian = [2*p*q^2*v
                              9*q*v^2 ];
            
            % Hessians
            if nargout >= 3
                % Hessian for p_k
                stateHessians(:, :, 1) = [   0    2*q*v^2
                                          2*q*v^2 2*p*v^2];
                
                % Hessian for q_k
                stateHessians(:, :, 2) = [2 0
                                          0 0];
                
                % Hessian for p_k
                noiseHessians(:, :, 1) = 2*p*q^2;
                
                % Hessian for q_k
                noiseHessians(:, :, 2) = 18*q*v;
            end
        end
    end
end

If you're estimator only requires the Jacobians of hk(xk,vk), i.e., the EKF, you do not have to implement the second-order derivatives. Simply compute and return solely the Jacobian when executing derivative(). Of course, when switching for example to the EKF2, your system model will not work until you also compute the Hessians.