D-Optimal Design for Parameter Estimation in Discrete-Time Nonlinear Dynamic Systems

Abstract

An optimal input design method for parameter estimation in a discrete-time nonlinear system is presented in the paper to improve the observability and identification precision of model parameters. Determinant of the information matrix is used as the criterion function which is generally a nonconvex function about the input signals to be designed. To avoid the locally optimizing problem, a randomized design method is proposed by which a globally optimizing test plan other than input signals may be obtained. Then the randomized design can be approximated by a nonrandomized design about optimal inputs. An iterative algorithm integrated with dynamic programming is given and verified by a numerical example on experimental design for self-calibration tests of ISP system.

1. Introduction

Model parameters applied to computation or compensation in science and engineering, such as error model coefficients in INS (inertial navigation system), generally require much higher identification precision than in other applications. However, haphazard experiments not only lead to poor accuracy in parameter estimating, but also would make some parameters unobservable. A good experimental design can increase both the precision and the efficiency of a test [1] and then improve the precision of system identification or state estimation [2–5].

The field of system identification and filtering are relatively mature [6–10]; relevant experimental design methods have not made substantial advance yet. D-optimal design which allocates the experimental input variables by maximizing the determinant of information matrix of the system is recognized as the most effective method for an experimental design [11, 12]. For model parameter identification of a dynamic system, the D-optimal design problem has a similar mathematical expression as the optimal control problem, but cannot be solved by the Pontryagin’s maximum principle and dynamic programming method, due to particularity in the form of performance index [13].

On the other hand, even for a dynamic system with linear or low-order nonlinear model, the D-optimal design problems may involve global optimizing of nonconvex function and cannot achieve analytical solutions by traditional nonlinear programming methods. Although many numerical searching algorithms have been proposed to solve the nonconvex problem in global optimizing, such as genetic algorithm, simulated annealing algorithm, and so forth, most of them are either time-consuming or no guarantee of global optimization of searching results [14–16].

D-optimal design for randomized inputs is a convex optimization technique, in which the experimental variables are transformed to test plans. A test plan specifies different probability measures to each input variable in admissible set and one selects inputs for a particular trial of the experiment via randomization. The randomized design method is mainly used in regression design problems. Mehra introduces the method to optimal input design for parameter identification in a discrete-time MIMO linear system with process noise [11].

Morelli and Klein consider input design problem for LTI systems in aircraft flight tests, and the specific goals with test time optimization are achieved using principles of dynamic programming [17]. Neto et al. generalize the results to nonlinear dynamic systems and consider additive colored noises in measurement [18]. The cost function selected by Neto et al. is the trace of a dispersion matrix in which the autocorrelation matrix of the colored measurement noises is introduced. They solve the optimization problem by genetic algorithm.

Lintereur studies optimal trajectories for a 2-axis gimbaled test table by which errors in inertial systems caused by angular motion are calibrated [19]. The trace of covariance matrix computed by Kalman filtering is minimized using a conjugate gradient algorithm, but local minimum may be obtained.

In this paper, we propose a randomized design method for parameter estimation in discrete-time nonlinear dynamic systems with constraints on inputs. By this design method, the original nonconvex optimization problem can be solved by the convex optimization technique, and the global optimal maximum is guaranteed. An iterative algorithm is given and verified by a numerical example on experimental design for continuous tumbling self-calibration tests of ISP system.

2. Problem Statement

In this section we give a mathematical formulation of the D-optimal design problem, in which a time-varying MIMO nonlinear system with unknown model parameters is considered. For simplification in notation and deduction, the process noise is assumed to be zero here.

Consider the following nonlinear dynamic system: 𝑥_𝑘+1=𝑓𝑥_𝑘,𝑢_𝑘,𝑘,𝜃,𝑥₀=𝑥₀,𝑦_𝑘=ℎ𝑥_𝑘,𝑡_𝑘,𝜃+𝑣_𝑘,𝑘=0,1,…,𝑁, (2.1) where 𝑥_𝑘=𝑥(𝑡_𝑘) is a 𝑛×1 state vector, 𝑥₀ is a constant vector, 𝑢_𝑘=𝑢(𝑡_𝑘) is a 𝑞×1 input vector, 𝑦_𝑘=𝑦(𝑡_𝑘) is a 𝑝×1 sampling output vector at moment k, and 𝑣_𝑘=𝑣(𝑡_𝑘) is a 𝑝×1 measurement noise vector. 𝑣_𝑘 is the Gaussian white noise sequences with 𝐸(𝑣_𝑘)=0 and 𝐸(𝑣_𝑘𝑣^𝑇_𝜏)=𝑅_𝑘𝛿_𝑘,𝜏, 𝑁 is the number of output samples observed and is fixed. Since 𝜃=[𝜃₁𝜃₂⋯𝜃_𝑚]^𝑇 denotes the 𝑚×1 vector of constant identifiable parameters, we estimate 𝜃 from the knowledge of {𝑦_𝑘,𝑢_𝑘−1,𝑘=1,…,𝑁} and give an unbiased efficient estimator 𝜃 with covariance 𝑀⁻¹, where 𝑀 is the Fisher information matrix. Therefore, the design problem is to select a series of inputs 𝑢_𝑘∈Ω_𝑢 such that a suitable criterion function corresponding to the objectives of the identification experiment is optimized.

The Fisher information matrix is defined as follows: 𝑀=𝐸_𝜃,𝑌𝜕log𝑝(𝑌,𝜃)()𝜕𝜃𝜕log𝑝𝑌,𝜃𝜕𝜃^T, (2.2) where 𝑌 denotes the set of observations {𝑦_𝑘,𝑘=1,…,𝑁} and the expectation in (2.2) is taken over by the sample space Ω_𝑌 of observations and the parameter space Ω_𝜃 of 𝜃.

Using conditional expectations, 𝑀 may be evaluated in two steps, first by computing 𝑀^′(𝜃)=𝐸_𝑌∣𝜃{•} and then 𝑀=𝐸_𝜃𝑀^′(𝜃). The second step is generally more tedious and an a priori distribution 𝑝(𝜃) should be known exactly. Here a Taylor-series approximation is used to simplify the computation: 𝑀^′_𝑖,𝑗(𝜃)=𝑀^′_𝑖,𝑗𝜃₀+𝜕𝑀^′_𝑖,𝑗𝜕𝜃_𝜃0𝜃−𝜃₀+12𝜃−𝜃₀^𝑇𝜕²𝑀^′_𝑖,𝑗𝜕𝜃²_𝜃0𝜃−𝜃₀+⋯, (2.3) where 𝜃₀ is the a priori mean of 𝜃,𝑖,𝑗=1,…,𝑚.

Retaining terms up to second order, 𝑀_𝑖,𝑗=𝑀^′_𝑖,𝑗𝜃₀+12𝜕tr²𝑀^′_𝑖,𝑗𝜕𝜃²_𝜃0𝑃₀, (2.4) where 𝑃₀ is the a priori covariance of 𝜃.

The second term is typically small compared to the first term either because 𝑃₀ is small or 𝑀^′(𝜃) is insensitive to 𝜃.

The conditional likelihood function 𝐿(𝜃)=log𝑝(𝑌∣𝜃) for system (2.1) is given as follows: 𝐿(𝜃)=−𝑁𝑝2()−1log2𝜋2^𝑁_𝑘=1𝑣_𝑘^𝑇𝑅_𝑘⁻¹𝑣_𝑘+log𝑅_𝑘. (2.5)

The matrix 𝑀^′(𝜃) has elements 𝑀^′_𝑖,𝑗(𝜃)=𝐸_𝑌𝜕𝐿(𝜃)𝜕𝜃_𝑖𝜕𝐿(𝜃)𝜕𝜃_𝑗=^𝑁_𝐾=1𝜕ℎ𝑥_𝑘,𝑡_𝑘,𝜃𝜕𝜃_𝑖^𝑇𝑅_𝑘⁻¹𝜕ℎ𝑥_𝑘,𝑡_𝑘,𝜃𝜕𝜃_𝑗. (2.6)

The sensitivity function is 𝜕ℎ𝑥_𝑘,𝑡_𝑘,𝜃𝜕𝜃_𝑖=𝜕ℎ𝜕𝑥^𝑇_𝑘𝑥_{𝜃𝑖,𝑘}+𝜕ℎ𝜕𝜃_𝑖, (2.7) where 𝑥_{𝜃𝑖,𝑘} is the partial derivative of 𝑥(𝑡) about 𝜃_𝑖 at moment k, that is, 𝑥_{𝜃𝑖,𝑘}=𝑥_𝜃𝑖(𝑡_𝑘) which meets the following equation: 𝑥_{𝜃𝑖,𝑘+1}=𝜕𝑓𝜕𝑥^𝑇_𝑘𝑥_{𝜃𝑖,𝑘}+𝜕𝑓𝜕𝜃_𝑖,𝑥_𝜃𝑖,0=0,𝑖=1,…,𝑚. (2.8)

Since 𝑥_{𝜃𝑖,𝑘} and 𝑥_𝑘 all depend on the elements of 𝑈, where 𝑈^T=[𝑢^𝑇₀⋯𝑢^𝑇_𝑁−1]∈Ω_𝑈 is the 𝑁𝑞-dimensional vector to be designed, we denote 𝑀^′(𝜃) as 𝑀^′(𝜃,𝑈).

From (2.6), it is easy to get 𝑀^′(𝜃,𝑈)=^𝑁_𝑘=1𝜕ℎ𝑥_𝑘,𝑡_𝑘,𝜃𝜕𝜃^𝑇𝑅_𝑘⁻¹𝜕ℎ𝑥_𝑘,𝑡_𝑘,𝜃.𝜕𝜃 (2.9)

Also from (2.4), the Fisher information matrix is generally 𝑀(𝑈)=𝐸_𝜃𝑀^′(𝜃,𝑈)≈𝑀^′𝜃₀,𝑈. (2.10)

There are many formulations of criterion function that measures the degree of observability about parameter 𝜃, such as tr(𝑀⁻¹(𝑈)) or |𝑀⁻¹(𝑈)|. A design which minimizes the scalar measure |𝑀⁻¹(𝑈)| or maximizes |𝑀(𝑈)| is called D-optimal, and it is equivalent to minimizing the volume of the uncertainty ellipsoid about parameter estimators. An important advantage of D-optimality is that it is invariant under scale changes in the parameters and linear transformations of the output.

Now, we choose |𝑀⁻¹(𝑈)| as the criterion function and formulate the D-optimal design problem as follows: min_{𝑈∈Ω𝑈}𝑀⁻¹(𝑈)s.t.𝑥_𝑘+1=𝑓𝑥_𝑘,𝑢_𝑘,𝑘,𝜃,𝑥₀=𝑥₀𝑥_{𝜃𝑖,𝑘+1}=𝜕𝑓𝜕𝑥^𝑇_𝑘𝑥_{𝜃𝑖,𝑘}+𝜕𝑓𝜕𝜃_𝑖,𝑥_𝜃𝑖,0=0,𝑖=1,…,𝑚. (2.11)

It should be pointed out that the problem in (2.11) cannot be solved by typical methods such as the Pontryagin’s maximum principle and dynamic programming method since |𝑀⁻¹(𝑈)| or tr(𝑀⁻¹(𝑈)) cannot be transformed to the index form in multistage decision process. In fact, there also exists great difficulty in getting a numerical solution with global optimization since |𝑀⁻¹(𝑈)| is not a convex function of 𝑈.

In the next section, we will present a randomized design method based on test plan theories.

3. D-Optimal Design Method

For a randomized input 𝑈∈Ω_𝑈 with probability measure 𝜉(𝑑𝑈)defined for all Borel sets and points of Ω_𝑈, the definition of the information matrix is 𝑀(𝜉)=_Ω𝑈𝑀(𝑈)⋅𝜉(𝑑𝑈), (3.1) where ∫_Ω𝑈𝜉(𝑑𝑈)=1.

If the probability measure is purely discrete, the information matrix is defined as follows: 𝑀(𝜉)=^𝑙_𝑖=1𝜉_𝑖𝑀𝑈_𝑖, (3.2) where 𝑙 is the number of spectrums, ∑^𝑙_𝑖=1𝜉_𝑖=1,  0≤𝜉_𝑖≤1.

𝑀(𝜉) is linear in 𝜉, so the criteria |𝑀⁻¹(𝜉)| or tr(𝑀⁻¹(𝜉)) are convex functions of 𝜉, and optimization with respect to 𝜉 gives globally optimizing design. However, we cannot find directly the optimal design 𝜉^∗ which minimizes |𝑀⁻¹(𝜉)|. Here, an iterative algorithm is proposed for searching 𝜉^∗ based on the following theorem [11].

Theorem 3.1. Let 𝜉^∗ be the optimal design then the following are equivalent:(i)𝜉^∗ maximizes |𝑀(𝜉)|,(ii)𝜉^∗ minimizes max_{𝑈∈Ω𝑈}tr(𝑀⁻¹(𝜉)𝑀(𝑈)),(iii)max_{𝑈∈Ω𝑈}tr(𝑀⁻¹(𝜉^∗)𝑀(𝑈))=tr(𝑀⁻¹(𝜉^∗)𝑀(𝜉^∗))=𝑚.
The D-optimal design 𝜉^∗ may be computed with the following algorithm.

Algorithm 3.2. Step 1. Start with any design 𝜉₀ such that 𝑀(𝜉₀) is nonsingular and let k=0.
Step 2. Compute 𝑀(𝜉_𝑘) and tr(𝑀⁻¹(𝜉_𝑘)𝑀(𝑈)) using (2.1), (2.6)–(2.10).
Step 3. Maximize tr(𝑀⁻¹(𝜉_𝑘)𝑀(𝑈)) over 𝑈∈Ω_𝑈 and get 𝑈_𝑘.
Step 4. If tr(𝑀⁻¹(𝜉_𝑘)𝑀(𝑈_𝑘))=𝑚, stop. Otherwise, let 𝜉_𝑘+1=(1−𝛼_𝑘)𝜉_𝑘+𝛼_𝑘𝜉(𝑈_𝑘),  0<𝛼_𝑘≤1 where 𝜉(𝑈_𝑘) is the design at the single point 𝑈_𝑘.
Choose 𝛼_𝑘 such that ∑^∞_𝑘=0𝛼_𝑘=∞, lim_𝑘→∞𝛼_𝑘=0, |𝑀(𝜉_𝑘+1)|≥|𝑀(𝜉_𝑘)|.
Step 5. Set 𝑘=𝑘+1 and go to step 2.
The convergence of the above algorithm to the global maximum is proved in the appendix.

Remarks. (1) The optimal design can be depicted by the following set: 𝜉₁,𝑈₁;𝜉₂,𝑈₂;…;𝜉_𝑙,𝑈_𝑙𝑚(),𝑙≤𝑚+12. (3.3) This can be used in a manner of randomized strategies when the experiment can be repeated. If the experiment is to be conducted only once, a nonrandomized design involving only one input 𝑈 should be preferable, that is, it assigns probability one to a particular 𝑈. Since the randomized design (3.3) has been derived, we can seek nonrandomized design 𝑈^∗ so that |𝑀(𝑈^∗)| approximates |∑^𝑙_𝑖=1𝜉_𝑖𝑀(𝑈_𝑖)|.
(2) Step  3 is most time-consuming computationally, and the criterion function tr(𝑀⁻¹(𝜉)𝑀(𝑈)) is generally not a convex function of 𝑈. Only if model (2.1) can be reduced to a linear discrete-time system, it would be a quadratic functional of 𝑈. Using (2.9) and (2.10), we get tr𝑀⁻¹(𝜉)𝑀(𝑈)≈tr𝑀⁻¹(𝜉)𝑀^′𝜃₀=,𝑈_𝑁𝑘=1tr𝑀⁻¹(𝜉)𝜕ℎ𝑥_𝑘,𝑡_𝑘,𝜃𝜕𝜃_𝜃0^𝑇𝑅_𝑘⁻¹𝜕ℎ𝑥_𝑘,𝑡_𝑘,𝜃𝜕𝜃_𝜃0. (3.4)

Unlike the computing of determinant, the operations with trace and sum of matrix can exchange order. By (3.4), the optimization problem can be solved by using maximum principle or dynamic programming methods since the above equation possesses the form of criterion function in multistage decision process.

Therefore, by using randomization and Theorem 3.1 the solution to a highly nonlinear and nonconvex optimization problem, that is, minimization of |𝑀⁻¹(𝑈)| is reduced to solving a relatively simpler optimization problem. This is mainly due to the fact that randomization produces convexity.

4. Simulation

In this section, we present a numerical example to verify the effectiveness of the proposed design method. The experiment to be designed is the continuous tumbling self-calibration test of an ISP system. Choose the determinant of information matrix as observability index and select the currents of gyro torquers or the command angular speed to the ISP as the experimental input variables, then the idea of D-optimal design can be applied to program the rotational trajectories of platform which represent the attitude and angular speed of platform at each moment.

First, we give the model equations of accelerometers and gyroscopes in the ISP system. The output equations of accelerometers that are also the observation equations are: 𝐲=0𝐴_𝑧−𝐴_𝑦−𝐴_𝑧0𝐴_𝑥𝐴_𝑦−𝐴_𝑥0𝜓_𝑥𝜓_𝑦𝜓_𝑧+𝑘_0𝑥+𝑘_1𝑥𝐴_𝑥𝑘_0𝑦+𝑘_1𝑦𝐴_𝑦−𝛼_𝑧𝐴_𝑥𝑘_0𝑧+𝑘_1𝑧𝐴_𝑧+𝛼_𝑦𝐴_𝑥−𝛼_𝑥𝐴_𝑦+𝜀, (4.1) where 𝛼_𝑥,𝛼_𝑦,𝛼_𝑧 represent the misalignment angles of accelerometers’ input axes with respect to platform frame.

𝑘_0𝑥,𝑘_0𝑦,𝑘_0𝑧 and 𝑘_1𝑥,𝑘_1𝑦,𝑘_1𝑧 represent the bias and scale factors of accelerometers.

𝑦 represents outputs of accelerometers; 𝜀 denotes the observation noises in outputs with zero mean and covariance matrix 𝛿²_𝑦𝐼_3×3 at each moment.

𝐴_𝑥, 𝐴_𝑦, 𝐴_𝑧 are projections of gravitational acceleration on the ideal platform frame which is defined based on 𝑥 accelerometer’s input axis and initially aligned to north, west, vertical direction, unit 𝑔₀, where 𝑔₀ is the local gravitational acceleration value: 𝐴_𝑥𝐴_𝑦𝐴_𝑧=(𝑎)(𝑐)(𝑎)(𝑏)(𝑐)(𝑎)(𝑐)(𝑎)(𝑏)(𝑐)(𝑎)(𝑏)sinsin−cossincossincos+cossinsincoscos, (4.2) where 𝑎,𝑏,𝑐 represent three ideal angular positions of platform gimbals from outer to inner, respectively, and meet the following differential equations: 𝑎̇𝑏̇𝑐̇=(𝑎)(𝑏)(𝑏)(𝑎)(𝑏)(𝑎)(𝑎)(𝑎)(𝑏)(𝑎)(𝑏)−costancossintan𝜔sin0coscossec0−sinsec^𝑒𝑐0𝜔_𝑒𝑠−(𝑐)(𝑏)(𝑐)(b)0(𝑐)(𝑐)0(𝑐)(𝑏)(𝑐)(𝑏)cossec−sinsecsincos−costansintan−1𝑡𝑔_𝑥𝑡𝑔_𝑦𝑡𝑔_𝑧, (4.3) where 𝜔_𝑒𝑐,𝜔_𝑒𝑠 represent north and vertical components of rotational speed of the earth 𝜔_𝑒, respectively.

𝜓_𝑥,𝜓_𝑦,𝜓_𝑧 in (4.1) represent the attitude errors between the practical platform frame and ideal one: 𝜓̇_𝑥𝜓̇_𝑦𝜓̇_𝑧=0𝑡𝑔_𝑧−𝑡𝑔_𝑦−𝑡𝑔_𝑧0𝑡𝑔_𝑥𝑡𝑔_𝑦−𝑡𝑔_𝑥0𝜓_𝑥𝜓_𝑦𝜓_𝑧+𝑑_0𝑥+𝑑_1𝑥𝑡𝑔_𝑥−𝛾_𝑥𝑧𝑡𝑔_𝑦+𝛾_𝑥𝑦𝑡𝑔_𝑧𝑑_0𝑦+𝑑_1𝑦𝑡𝑔_𝑦+𝛾_𝑦𝑧𝑡𝑔_𝑥−𝛾_𝑦𝑥𝑡𝑔_𝑧𝑑_0𝑧+𝑑_1𝑧𝑡𝑔_𝑧−𝛾_𝑧𝑦𝑡𝑔_𝑥+𝛾_𝑧𝑥𝑡𝑔_𝑦, (4.4) where 𝛾_𝑥𝑦,𝛾_𝑧𝑦,𝛾_𝑧𝑥,𝛾_𝑦𝑥,𝛾_𝑦𝑧,𝛾_𝑥𝑧 represent the misalignment angles of gyroscopes’ input axes with respect to platform frame.

𝑑_0𝑥,𝑑_0𝑦,𝑑_0𝑧 and 𝑑_1𝑥,𝑑_1𝑦,𝑑_1𝑧 represent the fixed drifts and scale factor errors of gyro torquers.

𝑢(𝑡)=[^{𝑡𝑔𝑥𝑡𝑔𝑦𝑡𝑔𝑧}]^𝑇 represent equivalent command angular speed to the ISP in an ideal platform frame and are the input variables to be designed.

Choose 𝜓_𝑥,𝜓_𝑦,𝜓_𝑧 and 𝑎,𝑏,𝑐 as the state vector 𝑥(𝑡); note that only the state of (4.4) is related to partial elements of unknown parameter vector 𝜃, where 𝜃^𝑇=𝜃^𝑇_𝑎𝜃^𝑇_𝑔,𝜃_𝑎=𝛼_𝑥,𝛼_𝑦,𝛼_𝑧,𝑘_0𝑥,𝑘_0𝑦,𝑘_0𝑧,𝑘_1𝑥,𝑘_1𝑦,𝑘_1𝑧^𝑇,𝜃_𝑔=𝛾_𝑥𝑦,𝛾_𝑧𝑦,𝛾_𝑧𝑥,𝛾_𝑦𝑥,𝛾_𝑦𝑧,𝛾_𝑥𝑧,𝑑_0𝑥,𝑑_0𝑦,𝑑_0𝑧,𝑑_1𝑥,𝑑_1𝑦,𝑑_1𝑧^T. (4.5)

Therefore, the equations for 𝑥_𝜃𝑖(𝑡) can be derived by making partial derivative of (4.4) about 𝜃_𝑔.

Rewrite (4.1), (4.3), and (4.4) for abbreviation as follows: 𝑥̇₁(𝑡)=𝐴₁(𝑢)𝑥₁(𝑡)+𝐵₁(𝑢)𝜃_𝑔,𝑥₁𝑡₀=𝑥₁₀,𝑥̇₂(𝑡)=𝑓𝑥₂,𝑢,𝑥₂𝑡₀𝑦=0,_𝑘=𝐻₁𝑥_2,𝑘𝑥_1,𝑘+𝐻₂𝑥_2,𝑘𝜃_𝑎+𝜀_𝑘,𝑘=0,1,…,𝑁, (4.6) where 𝑥₁(𝑡)=[𝜓_𝑥,𝜓_𝑦,𝜓_𝑧]^𝑇, 𝑥₂(𝑡)=[𝑎,𝑏,𝑐]^𝑇.

Then the equations for 𝑥_𝜃𝑖(𝑡) are as follows: 𝑥̇_𝜃𝑖(𝑡)=𝐴₁(𝑢)𝑥_𝜃𝑖(𝑡)+𝐵_1,𝑖(𝑢),𝑥_𝜃𝑖𝑡₀=0,𝑖=1,…,𝑠, (4.7) where 𝐵_1,𝑖(𝑢) denotes the 𝑖th column of matrix 𝐵₁(𝑢), and 𝑠 is the dimension of 𝜃_𝑔. 𝑀^′(𝜃,𝑈)=𝛿_𝑦^𝑁−2_𝑘=1𝐻^𝑇₂𝑥_2,𝑘𝐻₂𝑥_2,𝑘𝐻^𝑇₂𝑥_2,𝑘𝐻₁𝑥_2,𝑘𝑋_𝜃𝑋^𝑇_𝜃𝐻^𝑇₁𝑥_2,𝑘𝐻₂𝑥_2,𝑘𝑋^𝑇_𝜃𝐻^𝑇₁𝑥_2,𝑘𝐻₁𝑥_2,𝑘𝑋_𝜃 (4.8) Where 𝑋_𝜃=[𝑥_𝜃1,𝑘,…,𝑥_{𝜃𝑠,𝑘}].

It shows that 𝑀^′(𝜃,𝑈) is insensitive to 𝜃 since in (4.6) and (4.7), 𝑥_2,𝑘 as well as 𝑥_{𝜃𝑖,𝑘} is independent of the unknown parameter 𝜃.

Thus, in this problem 𝑀(𝑈)=𝑀^′(𝜃,𝑈),𝑀(𝜉)=^𝑙_𝑖=1𝜉_𝑖𝑀𝑈_𝑖 (4.9) for any 𝜃 and state equation for 𝑥₁(𝑡) can be deleted from the constraint conditions in (2.11).

Now we reformulate the design problem for self-calibration tests of ISP system as follows: max_{𝑈∈Ω𝑈}𝑀(𝜉𝑥̇)s.t.₂(𝑡)=𝑓𝑥₂,𝑢,𝑥₂𝑡₀𝑥̇=0,_𝜃𝑖(𝑡)=𝐴₁(𝑢)𝑥_𝜃𝑖(𝑡)+𝐵_1,𝑖(𝑢),𝑥_𝜃𝑖𝑡₀=0,𝑖=1,…,𝑠, (4.10) where Ω_𝑈={𝑈∶𝑢⁻≤𝑢_𝑘≤𝑢⁺,𝑘=0,…,𝑁−1} is the amplitude constraint on the precise command angular speed.

Since three state equations should be added to the constraint equations in (4.10) if one parameter in 𝜃_𝑔 is to be estimated, the dynamic programming algorithm in Step 3 will be very time-consuming. Therefore, only 9 parameters in 𝜃_𝑎 are considered in this example.

The initial design 𝜉₀ is chosen from a single-point design which maximizes |𝑀(𝑈)| via a rough search by confining each element of the input signals 𝑢(𝑡) to a four-segment square wave varying between −10 𝜔_𝑒, 0, and 10𝜔_𝑒. Figure 1 shows the convergence of tr(𝑀⁻¹(𝜉_𝑘)𝑀(𝑈_𝑘)) after eleven times of iterations, and |𝑀(𝜉_𝑘+1)|/|𝑀(𝜉_𝑘)| tends to one which means |𝑀(𝜉_𝑘)| converges to a local maximum |𝑀(𝜉^∗)|. Since |𝑀(𝜉)| is a concave function about 𝜉, 𝜉^∗ is global optimizing as well.

Figure 1 

Convergence of the proposed algorithm.

The values of log|𝑀(𝑈_𝑖)| at each supporting point 𝑈_𝑖 and measure 𝜉_𝑖 are listed in Table 1 with 𝛿²_𝑦=1. In addition, log|𝑀(𝜉^∗)|=log|∑_𝑖𝑀(𝑈_𝑖)𝜉_𝑖|=5.5 and log|𝑀(𝜉(𝑈′))|=log|𝑀(∑_𝑖𝑈_𝑖𝜉_𝑖)|=1.7. It shows that the randomized design 𝜉^∗ has better performance than any single-point design 𝜉(𝑈_𝑖), whereas the linear combination 𝑈′=∑_𝑖𝑈_𝑖𝜉_𝑖, although it is also an admissible control input in Ω_𝑈, shows poor performance and cannot be used as a nonrandomized approximation to the optimal test plan 𝜉^∗ here.

Table 1 

Design results by the propose algorithm.

In Table 1 the 4th supporting point 𝑈₄ has the maximum log|𝑀(𝑈_𝑖)| which approximates mostly to log|𝑀(𝜉^∗)|. So single-point design 𝜉(𝑈₄) is chosen as the nonrandomized approximation to 𝜉^∗. Comparison of estimation errors between tests with 𝑈₄ and 𝑈₁ (the initial rough design) is in Table 2, which shows some improvement in budgets of estimation precision. The control input curves of 𝑈₄ are plotted in Figure 2 with unit 𝜔_𝑒. The whole test time is three hours, and the sampling time is 22.5 minutes. None of 𝑡𝑔_𝑥, 𝑡𝑔_𝑦, and 𝑡𝑔_𝑧 in 𝑈₄ is the bang-bang type mainly because a constraint on angle 𝑏(−𝜋/6≤𝑏≤𝜋/6) is also assumed. The projection of gravitational acceleration in the ideal platform frame is plotted in Figure 3.

Table 2 

Standard errors of parameter estimators (𝛿_𝑦=1).

Figure 2 

Curves of the supporting point 𝑈4.

Figure 3 

Acceleration projection in the ideal platform frame.

The profiles of 𝐴_𝑥, 𝐴_𝑦, and 𝐴_𝑧 show that input axis of each accelerometer completes nearly a whole tumble in gravitational field which can provide sufficient stimulations to the error terms of accelerometers in practical testing.

5. Conclusion

A D-optimal design method for parameter estimation in nonlinear dynamic systems is presented based on test plan design theories. The corresponding iterative algorithm is proposed with a dynamic programming algorithm imbedded. The proof for the convergence of the algorithm is given as well. Simulation results on an optimal trajectory design problem in self-calibration test of ISP system demonstrate the effectiveness of the proposed algorithm.

Appendix

Proof of the Convergence of Algorithm 3.2

Proof. In the algorithm, 𝛼_𝑘 is chosen such that 𝑀𝜉₀≤𝑀𝜉₁≤⋯≤𝑀𝜉^∗. (A.1)
Since any bounded monotone nondecreasing sequence converges, the sequence |𝑀(𝜉₀)|,|𝑀(𝜉₁)|,…,|𝑀(𝜉_𝑘)| converges to some limit |𝑀(𝜉̂)|.
To prove |𝑀(𝜉̂)|=|𝑀(𝜉^∗)|, we assume the contrary |𝑀(𝜉̂)|<|𝑀(𝜉^∗)|. Then, by Theorem 3.1, for any 𝑘 there is a constant 𝜂 such that tr𝑀⁻¹𝜉_𝑘𝑀𝑈_𝑘−𝑚>𝜂>0. (A.2)
It follows that 𝜕𝜕𝛼_𝑘log𝑀𝜉_{𝑘+1𝛼𝑘=0}=tr𝑀⁻¹𝜉_𝑘−𝑀𝜉_𝑘+𝑀𝑈_𝑘>𝜂. (A.3)
For the smoothness of the function log|𝑀(𝜉_𝑘+1)| about 𝛼_𝑘 and by the assumption lim_𝑘→∞𝛼_𝑘=0, there is a positive integer 𝑠 such that for any 𝑘≥𝑠, (𝜕/𝜕𝛼_𝑘)log|𝑀(𝜉_𝑘+1)|≥𝜂.
Now integrating both sides of the above inequality over 𝛼_𝑘 from 0 to 𝛼_𝑘, one obtains 𝑀𝜉_𝑘+1𝑀𝜉_𝑘≥exp𝜂𝛼_𝑘. (A.4)
On the other hand, in view of the convergence of the sequence |𝑀(𝜉₀)|,…, for any small positive number 𝛾, there is a positive integer 𝑛 such that for any integer 𝑝≥𝑞≥𝑛, the following inequality holds: 𝑀𝜉_𝑝−𝑀𝜉_𝑞≤𝛾. (A.5)
Let 𝑡=max(𝑛,𝑠), then for any integer 𝑝≥𝑞≥𝑡, 𝛾≥𝑀𝜉_𝑝−𝑀𝜉_𝑞≥exp𝜂^𝑝−1_𝑘=𝑞𝛼_𝑘−1⋅𝑀𝜉_𝑞. (A.6)
This means ^𝑝−1_𝑘=𝑞𝛼_𝑘≤1𝜂log𝛾+𝑀𝜉_𝑞−log𝑀𝜉_𝑞, (A.7) which contradicts the assumption ∑^∞_𝑘=0𝛼_𝑘=∞.
Therefore, |𝑀(𝜉̂)|=|𝑀(𝜉^∗)|, the global maximum is obtained by the algorithm.

Acknowledgment

This work was supported by the Fundamental Research Funds for the Central Universities (Grant no. HIT. NSRIF. 2013037).