Adjustment of initial state vector (least squares)
[X0, info] = CL_adjustState(t0, X0, t, Zref, fun [[, opts]]) [X0, info] = CL_adjustState(t0, X0, t, Zref, list(fun, p1, p2 ... pn) [[, opts]])
Adjust the initial state vector so as to minimize measurements residuals.
The initial state vector is X0 (at time t0). The function fun propagates the state from time t0 to time t and computes he measurements vector Z(t) at time t. The quantity that is minimized in the sum of w_tk^2 * |Z_tk - Zref_tk|^2, where Z_tk and Zref_tk are respectively the kth components of Z(t) and Zref(t), and w_tk is the associated weight. Zref(t) represents the "real" measurements vector at time t.
The function fun must have the following interface:
[Z [, jac]] = fun(t0, X0, t [, p1, p2, ... pn])
With:
- t0: Initial time (1x1)
- X0: Initial state vector (at time t0) (Px1)
- t: Computation times (1xN)
- p1, p2, ... pn: Constant parameters of any type passed to the function
- Z: Computed measurements at times t (KxN). Undefined measurements components should be set to %nan, in order not to be considered in the computation.
- jac: Jacobian = dZ/dX0 (KxPxN). It may or not be computed. If it is not, it is safer to set it to [].
The structure opts contains adjustment options:
- indSf; (integer, size=1xQ) indices of solved-for parameters in state vector. Default is 1:P
- useJac: (boolean) Indicator of whether the Jacobian d(Z)/d(X0) should be used. %t => Jacobian is used, %f => Jacobian is not used; it is then evaluated automatically using finite difference. Default is %f
- svScale: Scale factors on the state vector components for a better matrix conditioning (size = Px1). Default is ones(P,1)
- eps: coefficient for finite difference calculation. State vector increments are equal to eps * svScale. Default is eps = 1.e-8
- derivOrder: Order for the computation of the Jacobian if computed using finite difference. Possible values are 1, 2 or 4. Default value is 1.
- wm: measurements weights for each type of measurement, that is each row of the measurement vector (Kx1). In theory, the weight should be equal to 1/sigmak where sigmak is the standard deviation of the error on measurement k (k in 1 .. K). Default value is ones(K,1)
- tau: time constant for additional exponential weighting. Default is %inf (no exponential weighting)
- residMargin: Criterion for iterations control. Iterations stop when the relative variation of residuals (RMS) between 2 iterations is less than this value. Default is 1.e-3 (0.1%).
- verbose: Verbose level. 0 means no information is printed on the screen, 1 means information is printed on the screen. Default is 0.
- meth: Method used to compute state vector adjustments at each iteration. Can be "lsq" (Scilab "lsq" function is used) or "pinv" (Scilab function "pinv" is used). Default is "pinv".
- nIterMax: Maximum number of iterations. Default value is 30.
- nTrials: Number of successive cases for which residuals don't decrease (including margin). Default value is 3.
info is a structure that contains additional information:
- niter: number of iterations before the algorithm stops
- residuals_rms: root mean square of residuals at the last iteration
- res constains detailed information for the last iteration, in particular:
- X0: Initial state vector before adjustment
- resid_rms: Residuals RMS corresponding to X0
- X0_adj: Initial state vector after adjustment
The algorithm stops if the following conditions are met:
- the residuals (RMS) are greater than the minimum value obtained 'nTrials' times in a row (margin included), or
- the maximum number of iterations has been reached, or
- the residuals are exactly 0.
Initial time (1x1)
Initial/adjusted state vector (Px1)
Times at which measurements are given (1xN)
Reference measurements at times t (KxN)
Function that computes the measurements: function or list(function [, p1, ...]). See above for details.
(structure, optional) Computation options. See above for details.
(structure) Additional information on the results. See above for details.
CNES - DCT/SB
// ------------------------------ // Case 1: Jacobian is not computed // ------------------------------ // Model function [z]=F1(t0, x0, t, p) z = x0(1) + x0(2)^2 * (t - t0) + p; endfunction // Initial and measurements times t0 = 0; t = t0 + linspace(0, 2, 100); p1 = 0; // Exact solution x0_ref = [1; 2]; z_ref = F1(t0, x0_ref, t, p1); // Initial state vector estimate x0_est = x0_ref + [-0.2; 0.1]; // Options: default values // "useJac" explicitly set although %f is the default opts = struct(); opts.useJac = %f; x0 = CL_adjustState(t0, x0_est, t, z_ref, list(F1, p1), opts); x0 - x0_ref // ------------------------------ // Case 2: Jacobian is computed // and used // ------------------------------ function [z, jac]=F2(t0, x0, t, p) z = F1(t0, x0, t, p); jac = zeros(1, 2, length(t)); jac(1, 1, :) = 1; jac(1, 2, :) = 2 * x0(2) * (t - t0); endfunction // Options: change value of "useJac" opts.useJac = %t; x0 = CL_adjustState(t0, x0_est, t, z_ref, list(F2, p1), opts); x0 - x0_ref | ![]() | ![]() |