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

Authors: Yu Liu; Hamid Reza Karimi; Zhiwei Yu
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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:(2.1)xk+1=f(xk,uk,k,θ),x0=x-0,yk=h(xk,tk,θ)+vk,k=0,1,…,N,
where xk=x(tk) is a n×1 state vector, x-0 is a constant vector, uk=u(tk) is a q×1 input vector, yk=y(tk) is a p×1 sampling output vector at moment k, and vk=v(tk) is a  p×1 measurement noise vector. vk is the Gaussian white noise sequences with E(vk)=0 and E(vkvτT)=Rkδk,τ, N is the number of output samples observed and is fixed. Since  θ=[θ1θ2⋯θm]T denotes the m×1 vector of constant identifiable parameters, we estimate θ from the knowledge of {yk,uk-1,k=1,…,N} and give an unbiased efficient estimator θ^  with covariance M-1, where M is the Fisher information matrix. Therefore, the design problem is to select a series of inputs uk∈Ωu such that a suitable criterion function corresponding to the objectives of the identification experiment is optimized.The Fisher information matrix is defined as follows:(2.2)M=Eθ,Y[(∂logp(Y,θ)∂θ)(∂logp(Y,θ)∂θ)T],
where Y denotes the set of observations {yk,k=1,…,N} and the expectation in (2.2) is taken over by the sample space ΩY of observations and the parameter space Ωθ of θ.Using conditional expectations,M may be evaluated in two steps, first by computing M′(θ)=EY∣θ{•} and then M=EθM′(θ). The second step is generally more tedious and an a priori distribution p(θ) should be known exactly. Here a Taylor-series approximation is used to simplify the computation:
(2.3)M′i,j(θ)=M′i,j(θ0)+∂Mi,j′∂θ|θ0(θ-θ0)+12(θ-θ0)T∂2Mi,j′∂θ2|θ0(θ-θ0)+⋯,
where θ0 is the a priori mean of  θ,i,j=1,…,m.Retaining terms up to second order,(2.4)Mi,j=Mi,j′(θ0)+12tr[∂2M′i,j∂θ2|θ0P0],
where P0 is the a priori covariance of θ.The second term is typically small compared to the first term either becauseP0 is small or M′(θ) is insensitive to θ.The conditional likelihood functionL(θ)=logp(Y∣θ) for system (2.1) is given as follows:
(2.5)L(θ)=-Np2log(2π)-12∑k=1N{(vk)TRk-1vk+log|Rk|}.The matrixM′(θ) has elements
(2.6)Mi,j′(θ)=EY(∂L(θ)∂θi∂L(θ)∂θj)=∑K=1N{(∂h(xk,tk,θ)∂θi)TRk-1∂h(xk,tk,θ)∂θj}.The sensitivity function is(2.7)∂h(xk,tk,θ)∂θi=∂h∂xkTxθi,k+∂h∂θi,
where xθi,k is the partial derivative of x(t) about θi at moment k, that is, xθi,k=xθi(tk) which meets the following equation:
(2.8)xθi,k+1=∂f∂xkTxθi,k+∂f∂θi,xθi,0=0,i=1,…,m.Sincexθi,k and xk all depend on the elements of U, where UT=[u0T⋯uN-1T]∈ΩU is the Nq-dimensional vector to be designed, we denote M′(θ) as M′(θ,U).From (2.6), it is easy to get
(2.9)M′(θ,U)=∑k=1N{(∂h(xk,tk,θ)∂θ)TRk-1∂h(xk,tk,θ)∂θ}.Also from (2.4), the Fisher information matrix is generally
(2.10)M(U)=EθM′(θ,U)≈M′(θ0,U).There are many formulations of criterion function that measures the degree of observability about parameterθ, such as tr(M-1(U)) or |M-1(U)|. A design which minimizes the scalar measure |M-1(U)| or maximizes |M(U)| 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|M-1(U)| as the criterion function and formulate the D-optimal design problem as follows:
(2.11)minU∈ΩU|M-1(U)|s.t.xk+1=f(xk,uk,k,θ),x0=x-0s.t.xθi,k+1=∂f∂xkTxθi,k+∂f∂θi,xθi,0=0,i=1,…,m.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 |M-1(U)| or tr(M-1(U)) 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 |M-1(U)| is not a convex function of  U.In the next section, we will present a randomized design method based on test plan theories.


## 3. D-Optimal Design Method

For a randomized inputU∈ΩU with probability measure ξ(dU)defined for all Borel sets and points of ΩU, the definition of the information matrix is
(3.1)M(ξ)=∫ΩUM(U)·ξ(dU),
where ∫ΩUξ(dU)=1.If the probability measure is purely discrete, the information matrix is defined as follows:(3.2)M(ξ)=∑i=1lξiM(Ui),
where l is the number of spectrums,  ∑i=1lξi=1,  0≤ξi≤1.M
(
ξ
) is linear in ξ, so the criteria |M-1(ξ)| or tr(M-1(ξ)) are convex functions of ξ, and optimization with respect to ξ gives globally optimizing design. However, we cannot find directly the optimal design ξ* which minimizes |M-1(ξ)|. 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 |M(ξ)|,(ii)
ξ* minimizes maxU∈ΩUtr(M-1(ξ)M(U)),(iii)
maxU∈ΩUtr(M-1(ξ*)M(U))=tr(M-1(ξ*)M(ξ*))=m.
The D-optimal designξ* may be computed with the following algorithm.Algorithm 3.2.
Step 1. Start with any design ξ0 such that M(ξ0) is nonsingular and let k=0.
Step 2. Compute M(ξk) and  tr(M-1(ξk)M(U)) using (2.1), (2.6)–(2.10).
Step 3. Maximize tr(M-1(ξk)M(U)) over U∈ΩU  and  get Uk.
Step 4. If tr(M-1(ξk)M(Uk))=m, stop. Otherwise, let ξk+1=(1-αk)ξk+αkξ(Uk),  0<αk≤1 where ξ(Uk) is the design at the single point Uk.
Chooseαk such that ∑k=0∞αk=∞,  limk→∞αk=0, |M(ξk+1)|≥|M(ξk)|.
Step 5. Set k=k+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:
(3.3){ξ1,U1;ξ2,U2;…;ξl,Ul},l≤m(m+1)2.
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 U should be preferable, that is, it assigns probability one to a particular U. Since the randomized design (3.3) has been derived, we can seek nonrandomized design U* so that |M(U*)| approximates |∑i=1lξiM(Ui)|.
(
2
) Step  3 is most time-consuming computationally, and the criterion function tr(M-1(ξ)M(U)) is generally not a convex function of  U. Only if model (2.1) can be reduced to a linear discrete-time system, it would be a quadratic functional of U. Using (2.9) and (2.10), we get
(3.4)tr(M-1(ξ)M(U))≈tr(M-1(ξ)M′(θ0,U))=∑k=1Ntr{M-1(ξ)(∂h(xk,tk,θ)∂θ|θ0)TRk-1∂h(xk,tk,θ)∂θ|θ0}.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 Theorem3.1 the solution to a highly nonlinear and nonconvex optimization problem, that is, minimization of |M-1(U)| 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:(4.1)y=[0Az-Ay-Az0AxAy-Ax0][ψxψyψz]+[k0x+k1xAxk0y+k1yAy-αzAxk0z+k1zAz+αyAx-αxAy]+ε,
where αx,αy,αz represent the misalignment angles of accelerometers’ input axes with respect to platform frame.k
0
x
,
k
0
y
,
k
0
z  and  k1x,k1y,k1z represent the bias and scale factors of accelerometers.y represents outputs of accelerometers; ε denotes the observation noises in outputs with zero mean and covariance matrix  δy2I3×3  at each moment.A
x, Ay, Az are projections of gravitational acceleration on the ideal platform frame which is defined based on x accelerometer’s input axis and initially aligned to north, west, vertical direction, unit g0, where g0 is the local gravitational acceleration value:
(4.2)[AxAyAz]=[sin(a)sin(c)-cos(a)sin(b)cos(c)sin(a)cos(c)+cos(a)sin(b)sin(c)cos(a)cos(b)],
where a,b,c represent three ideal angular positions of platform gimbals from outer to inner, respectively, and meet the following differential equations:
(4.3)[a˙b˙c˙]=[-cos(a)tan(b)cos(b)sin(a)tan(b)sin(a)0cos(a)cos(a)sec(b)0-sin(a)sec(b)][ωec0ωes]-[cos(c)sec(b)-sin(c)sec(b)0sin(c)cos(c)0-cos(c)tan(b)sin(c)tan(b)-1][tgxtgytgz],
where ωec,ωes represent north and vertical components of rotational speed of the earth ωe, respectively.ψ
x
,
ψ
y
,
ψ
z in (4.1) represent the attitude errors between the practical platform frame and ideal one:
(4.4)[ψ˙xψ˙yψ˙z]=[0tgz-tgy-tgz0tgxtgy-tgx0][ψxψyψz]+[d0x+d1xtgx-γxztgy+γxytgzd0y+d1ytgy+γyztgx-γyxtgzd0z+d1ztgz-γzytgx+γzxtgy],
where γxy,γzy,γzx,γyx,γyz,γxz represent the misalignment angles of gyroscopes’ input axes with respect to platform frame.d
0
x
,
d
0
y
,
d
0
z and d1x,d1y,d1z represent the fixed drifts and scale factor errors of gyro torquers.u
(
t
)
=
[
t
g
x
t
g
y
t
g
z
]
T  represent equivalent command angular speed to the ISP in an ideal platform frame and are the input variables to be designed.Chooseψx,ψy,ψz and a,b,c as the state vector x(t); note that only the state of (4.4) is related to partial elements of unknown parameter vector θ, where
(4.5)θT=[θaTθgT],θa=[αx,αy,αz,k0x,k0y,k0z,k1x,k1y,k1z]T,θg=[γxy,γzy,γzx,γyx,γyz,γxz,d0x,d0y,d0z,d1x,d1y,d1z]T.Therefore, the equations forxθi(t) can be derived by making partial derivative of (4.4) about θg.Rewrite (4.1), (4.3), and (4.4) for abbreviation as follows:
(4.6)x˙1(t)=A1(u)x1(t)+B1(u)θg,x1(t0)=x10,x˙2(t)=f(x2,u),x2(t0)=0,yk=H1(x2,k)x1,k+H2(x2,k)θa+εk,k=0,1,…,N,
where  x1(t)=[ψx,ψy,ψz]T, x2(t)=[a,b,c]T.Then the equations forxθi(t) are as follows:
(4.7)x˙θi(t)=A1(u)xθi(t)+B1,i(u),xθi(t0)=0,i=1,…,s,
where B1,i(u) denotes the ith column of matrix B1(u), and s is the dimension of  θg.
(4.8)M′(θ,U)=δy-2∑k=1N[H2T(x2,k)H2(x2,k)H2T(x2,k)H1(x2,k)XθXθTH1T(x2,k)H2(x2,k)XθTH1T(x2,k)H1(x2,k)Xθ]
Where  Xθ=[xθ1,k,…,xθs,k].It shows thatM′(θ,U) is insensitive to θ since in (4.6) and (4.7), x2,k as well as xθi,k is independent of the unknown parameter θ.Thus, in this problem(4.9)M(U)=M′(θ,U),M(ξ)=∑i=1lξiM(Ui)
for any θ and state equation for x1(t) 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:(4.10)maxU∈ΩU|M(ξ)|s.t.x˙2(t)=f(x2,u),x2(t0)=0,s.t.x˙θi(t)=A1(u)xθi(t)+B1,i(u),xθi(t0)=0,i=1,…,s,
where ΩU={U:u-≤uk≤u+,k=0,…,N-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 θg is to be estimated, the dynamic programming algorithm in Step 3 will be very time-consuming. Therefore, only 9 parameters in θa are considered in this example.The initial designξ0 is chosen from a single-point design which maximizes |M(U)| via a rough search by confining each element of the input signals u(t) to a four-segment square wave varying between −10 ωe, 0, and 10ωe. Figure 1 shows the convergence of tr(M-1(ξk)M(Uk)) after eleven times of iterations, and |M(ξk+1)|/|M(ξk)| tends to one which means |M(ξk)| converges to a local maximum |M(ξ*)|. Since |M(ξ)| is a concave function about ξ, ξ* is global optimizing as well.Figure 1
Convergence of the proposed algorithm.The values oflog|M(Ui)| at each supporting point Ui and measure ξi are listed in Table 1 with δy2=1. In addition, log|M(ξ*)|=log|∑iM(Ui)ξi|=5.5 and log|M(ξ(U'))|=log|M(∑iUiξi)|=1.7. It shows that the randomized design ξ* has better performance than any single-point design ξ(Ui), whereas the linear combination U'=∑iUiξi, although it is also an admissible control input in ΩU, 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.
i
1
2
3
4
log
|
M
(
U
i
)
|
4.7
4.5
4.7
5.2
ξ
i
0.1
0.2
0.2
0.5In Table1 the 4th supporting point U4 has the maximum log|M(Ui)| which approximates mostly to log|M(ξ*)|. So single-point design ξ(U4) is chosen as the nonrandomized approximation to ξ*. Comparison of estimation errors between tests with U4 and U1 (the initial rough design) is in Table 2, which shows some improvement in budgets of estimation precision. The control input curves of U4 are plotted in Figure 2 with unit ωe. The whole test time is three hours, and the sampling time is 22.5 minutes. None of  tgx, tgy, and tgz in U4 is the bang-bang type mainly because a constraint on angle b(-π/6≤b≤π/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 (δy=1).
α
x
α
y
α
z
k
0
x
k
0
y
k
0
z
k
1
x
k
1
y
k
1
z
U
1
0.67
0.77
0.76
0.4
0.43
0.43
0.7
0.67
0.63
U
4
0.65
0.57
0.53
0.36
0.36
0.38
0.53
0.56
0.66Figure 2
Curves of the supporting pointU4.Figure 3
Acceleration projection in the ideal platform frame.The profiles ofAx, Ay, and Az 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.



---

*Source: 296701-2012-09-20.xml*