ABSTRACT

In this paper, an efficient computational approach is proposed to solve the discrete time nonlinear stochastic optimal control problem. For this purpose, a linear quadratic regulator model, which is a linear dynamical system with the quadratic criterion cost function, is employed. In our approach, the model-based optimal control problem is reformulated into the input-output equations. In this way, the Hankel matrix and the observability matrix are constructed. Further, the sum squares of output error is defined. In these point of views, the least squares optimization problem is introduced, so as the differences between the real output and the model output could be calculated. Applying the first-order derivative to the sum squares of output error, the necessary condition is then derived. After some algebraic manipulations, the optimal control law is produced. By substituting this control policy into the input-output equations, the model output is updated iteratively. For illustration, an example of the direct current and alternating current converter problem is studied. As a result, the model output trajectory of the least squares solution is close to the real output with the smallest sum squares of output error. In conclusion, the efficiency and the accuracy of the approach proposed are highly presented.

Keywords:

Least Squares Solution, Stochastic Optimal Control, Linear Quadratic Regulator, Sum Squares of Output Error, Input-Output Equations

1. Introduction

Stochastic dynamical system is a practical system in modeling and simulating the real-world problems. The behavior of the fluctuation, which is caused by the effect of noise disturbance in the dynamical system to represent the real situation, rises to the attention of many researchers. See for examples, [1] - [8] . As such, optimization and control of the stochastic dynamical system are a challenging topic. In the past decades, the research outcomes, both for theoretical results and development of the algorithms, are well-defined [9] [10] [11] [12] [13] . In general, the real processes, which are modeled into the stochastic optimal control problem, are mostly nonlinear process. Their actual models in natural could be necessary unknown. Applying the linear model, including the Kalman filtering theory, to solve the stochastic optimal control problems, in particular, for nonlinear case, give a simplification methodology instead [14] [15] [16] [17] .

Recently, the integrated optimal control and parameter estimation (ICOPE) algorithm, which solves the linear model-based optimal control problem iteratively, is proposed in the literature [18] [19] . The purpose of this algorithm is to provide the optimal solution of the discrete time nonlinear stochastic optimal control problem with the different structure and parameters. In this algorithm, the adjusted parameters are added into the model used so as the differences between the real plant and the model used can be calculated, in turn, update the optimal solution of the model used during the computation procedure. The integration of system optimization and parameter estimation is the main feature of this algorithm, where the principle of model-reality differences and the Kalman filtering theory are incorporated. Once the convergence is achieved, the iterative solution approximates to the correct optimal solution of the original optimal control problem, in spite of the model-reality differences. However, increasing the accuracy of the model output to track the real output would be necessary required.

In this paper, an efficient matching scheme, which diminishes the adjusted parameters, is established deliberately. In our approach, the model-based optimal control problem is simplified from the discrete-time nonlinear stochastic optimal control problem. Then, this model-based optimal control problem is reformulated into the input-output equations. During the formulation of the input-output equations, the Hankel matrix and the observability matrix are constructed. These matrices capture the characteristic of the model used into the output measurement. By virtue of this, the least square optimization problem is introduced. From the validation of the first order necessary condition, the normal equation is resulted and the optimal control law is updated accordingly on the recursion formula. As a result of this, the sum squares of the output error could be minimized demonstratively. Hence, the efficiency of the algorithm proposed is highly presented.

The rest of the paper is organized as follows. In Section 2, the discrete time nonlinear stochastic optimal control problem and its simplified model-based optimal control problem are described. In Section 3, the system optimization with the least squares updating scheme is discussed. The computation procedure is summarized as the iterative algorithm. In Section 4, the current converter problem is illustrated and the result is demonstrated. Finally, some concluding remarks are made.

2. Problem Description

Consider a general class of dynamical system given by

$x\left(k+1\right)=f\left(x\left(k\right),u\left(k\right),k\right)+G\omega \left(k\right)$ (1a)

$y\left(k\right)=h\left(x\left(k\right),k\right)+\eta \left(k\right)$ (1b)

where $u\left(k\right)\in {\Re }^m,k=0,1,\cdots ,N-1,x\left(k\right)\in {\Re }^n,k=0,1,\cdots ,N$ , and

$y\left(k\right)\in {\Re }^p,k=0,1,\cdots ,N$ , are, respectively, the control sequence, the state sequence and the output sequence. The terms $\omega \left(k\right)\in {\Re }^q,k=0,1,\cdots ,N-1$ , and $\eta \left(k\right)\in {\Re }^p,k=0,1,\cdots ,N$ , are, respectively, process noise sequences and output noise sequences. Both of these noise sequences are the stationary Gaussian white noise sequences with zero mean and their covariance matrices are given by $Q_{\omega }$ and $R_{\eta }$ , respectively, where $Q_{\omega }$ is a $q\times q$ positive definite matrix and $R_{\eta }$ is a $p\times p$ positive definite matrix. In addition, $G$ is an $n\times q$ process noise coefficient matrix, $f:{\Re }^n\times {\Re }^m\times \Re \to {\Re }^n$ represents the plant dynamics and $h:{\Re }^n\times \Re \to {\Re }^p$ is the output measurement channel.

Here, the aim is to find the control sequence $u\left(k\right),k=0,1,\cdots ,N-1$ , such that the following cost function

$g_0\left(u\right)=E\left[\phi \left(x\left(N\right),N\right)+{\displaystyle \underset{k=0}{\overset{N-1}{\sum }}L\left(x\left(k\right),u\left(k\right),k\right)}\right]$ (2)

is minimized over the dynamical system in Equation (1), where $\phi :{\Re }^n\times \Re \to \Re$ is the terminal cost and $L:{\Re }^n\times {\Re }^m\times \Re \to \Re$ is the cost under summation. The cost function $g_0$ is the scalar function and $E[\cdot ]$ is the expectation operator. It is assumed that all functions in Equations (1) and (2) are continuously differentiable with respect to their respective arguments.

The initial state is

$x\left(0\right)=x_0$

where $x_0\in {\Re }^n$ is a random vector with mean and variance given, respectively, by

$E\left[x_0\right]={\bar{x}}_0$ and $E\left[\left(x_0-{\bar{x}}_0\right){\left(x_0-{\bar{x}}_0\right)}^{\text{T}}\right]=M_0$ .

Here, $M_0$ is an $n\times n$ positive definite matrix. It is assumed that the initial state, the process noise and the measurement noise are statistically independent.

This problem is regarded as the discrete-time nonlinear stochastic optimal control problem and is referred to as Problem (P). Notice that the structure of Problem (P) is complex, and the exact solution of Problem (P) is, in general, unable to be obtained. In view of these, Problem (P) is proposed to be solved via solving a simplified model-based optimal control problem iteratively. For this purpose, let this simplified model-based optimal control problem, which is referred to as Problem (M), be given below:

$\underset{u\left(k\right)}{\min }{g}_1\left(u\right)=\frac{1}{2}\bar{x}{\left(N\right)}^{\text{T}}S\left(N\right)\bar{x}\left(N\right)+\frac{1}{2}{\displaystyle \underset{k=0}{\overset{N-1}{\sum }}\left(\bar{x}{\left(k\right)}^{\text{T}}Q\bar{x}\left(k\right)+u{\left(k\right)}^{\text{T}}Ru\left(k\right)\right)}$

subject to (3)

$\bar{x}\left(k+1\right)=A\bar{x}\left(k\right)+Bu\left(k\right),$ $\bar{x}\left(0\right)={\bar{x}}_0$

$\bar{y}\left(k\right)=C\bar{x}\left(k\right)$

where $\bar{x}\left(k\right)\in {\Re }^n,k=0,1,\cdots ,N$ , and $\bar{y}\left(k\right)\in {\Re }^p,k=0,1,\cdots ,N$ , are, respectively, the expected state sequence and the expected output sequence. Here, $A$ is an $n\times n$ state transition matrix, $B$ is an $n\times m$ control coefficient matrix, $C$ is a $p\times n$ output coefficient matrix, whereas $S\left(N\right)$ and $Q$ are $n\times n$ positive semi-definite matrices, $R$ is a $m\times m$ positive definite matrix and $g_1$ is the scalar cost function.

Notice that only solving Problem (M) actually would not give the optimal solution of Problem (P). However, by establishing an efficient matching scheme, it is possible to approximate the true optimal solution of Problem (P), in spite of model-reality differences. This could be done iteratively.

3. System Optimization with Least Square Updating Scheme

Now, let us define the Hamiltonian function for Problem (M) as follows:

$H\left(k\right)=\frac{1}{2}\left(\bar{x}{\left(k\right)}^{\text{T}}Q\bar{x}\left(k\right)+u{\left(k\right)}^{\text{T}}Ru\left(k\right)\right)+p{\left(k+1\right)}^{\text{T}}\left(A\bar{x}\left(k\right)+Bu\left(k\right)\right)$ . (4)

Then, from Equation (3), the augmented cost function becomes

$\begin{array}{l}{g^{\prime }}_1\left(u\right)=\frac{1}{2}\bar{x}{\left(N\right)}^{\text{T}}S\left(N\right)\bar{x}\left(N\right)+p{\left(0\right)}^{\text{T}}\bar{x}\left(0\right)\\ -p{\left(N\right)}^{\text{T}}\bar{x}\left(N\right)+{\displaystyle \underset{k=0}{\overset{N-1}{\sum }}\left(H\left(k\right)-p{\left(k\right)}^{\text{T}}\bar{x}\left(k\right)\right)}\end{array}$ (5)

where $p\left(k\right)\in {\Re }^n$ is the appropriate multiplier to be determined later.

Applying the calculus of variation [20] [21] [22] to the augmented cost function in Equation (5), the following necessary conditions are resulted:

1) Stationary condition:

$\frac{\partial H\left(k\right)}{\partial u\left(k\right)}=Ru\left(k\right)+B^{\text{T}}p\left(k+1\right)=0$ (6)

2) Costate equation:

$\frac{\partial H\left(k\right)}{\partial \bar{x}\left(k\right)}=Q\bar{x}\left(k\right)+A^{\text{T}}p\left(k+1\right)=p\left(k\right)$ (7)

3) State equation:

$\frac{\partial H\left(k\right)}{\partial p\left(k+1\right)}=A\bar{x}\left(k\right)+Bu\left(k\right)=\bar{x}\left(k+1\right)$ (8)

with the boundary conditions $\bar{x}\left(0\right)={\bar{x}}_0$ and $p\left(N\right)=S\left(N\right)\bar{x}\left(N\right)$ .

3.1. State Feedback Optimal Control Law

The following theorem expresses the feedback control law based on the state that can be used in solving Problem (M).

Theorem 1. For the given Problem (M), the optimal control law is the feedback control law defined by

$u\left(k\right)=-K\left(k\right)\bar{x}\left(k\right)$ (9)

where

$K\left(k\right)={\left(B^{\text{T}}S\left(k+1\right)B+R\right)}^{-1}{B}^{\text{T}}S\left(k+1\right)A$ (10)

$S\left(k\right)=A^{\text{T}}S\left(k+1\right)\left(A-BK\left(k\right)\right)+Q$ (11)

with the boundary condition $S\left(N\right)$ given. Here, the feedback control law is a linear combination of the states. That is, the optimal control is linear state-vari- able feedback.

Proof. From Equation (6), the stationary condition can be rewritten as

$Ru\left(k\right)=-B^{\text{T}}p\left(k+1\right)$ . (12)

Applying the sweep method [20] [21] [22] ,

$p\left(k\right)=S\left(k\right)\bar{x}\left(k\right)$ , (13)

for $k=k+1$ in Equation (12) to yield

$Ru\left(k\right)=-B^{\text{T}}S\left(k+1\right)\bar{x}\left(k+1\right)$ . (14)

Taking Equation (8) into Equation (14), we have

$Ru\left(k\right)=-B^{\text{T}}S\left(k+1\right)\left(A\bar{x}\left(k\right)+Bu\left(k\right)\right)$ .

After some algebraic manipulations, the feedback control law (9) is obtained, where Equation (10) is satisfied.

Now, substituting Equation (13) for $k=k+1$ into Equation (7), the costate equation is written as

$p\left(k\right)=Q\bar{x}\left(k\right)+A^{\text{T}}S\left(k+1\right)\bar{x}\left(k+1\right)$ (15)

and considering the state Equation (8) in Equation (15), the costate equation becomes

$p\left(k\right)=Q\bar{x}\left(k\right)+A^{\text{T}}S\left(k+1\right)\left(A\bar{x}\left(k\right)+Bu\left(k\right)\right)$ . (16)

Hence, by applying the feedback control law (9) in Equation (16), and doing some algebraic manipulations, it could be seen that Equation (11) is satisfied after comparing the manipulation result to Equation (13). This completes the proof.

By substituting Equation (9) into Equation (8), the state equation becomes

$\bar{x}\left(k+1\right)=\left(A-BK\left(k\right)\right)\bar{x}\left(k\right)$ (17)

and the model output is measured from

$\bar{y}\left(k\right)=C\bar{x}\left(k\right)$ . (18)

In view of this, the calculation procedure for obtaining the feedback control law for Problem (M) is summarized below:

Algorithm 1: Feedback control algorithm

Data Given $A,B,C,Q,R,S\left(N\right),x_0,N$ .

Step 0 Calculate $K\left(k\right),k=0,1,\cdots ,N-1$ , and $S\left(k\right),k=0,1,\cdots ,N$ , from Equations (10) and (11), respectively.

Step 1 Solve Problem (M) that is defined by Equation (3) to obtain

$u\left(k\right),k=0,1,\cdots ,N-1$ , and $\bar{x}\left(k\right),\bar{y}\left(k\right),k=0,1,\cdots ,N$ , respectively, from Equations (9), (17) and (18).

Step 2 Evaluate the cost function $g_1$ from Equation (3).

Remarks:

1) Data A, B, C can be obtained by the linearization of the real plant $f$ and the output measurement $h$ from Problem (P).

2) In Step 0, the offline calculation is done for $K\left(k\right),k=0,1,\cdots ,N-1$ , and $S\left(k\right),k=0,1,\cdots ,N$ .

3) The solution procedure, which the dynamical system is solved in Step 1, and the cost function is evaluated in Step 2, is known as system optimization.

3.2. Least Squares Updating Scheme

Now, we define the output error $r:{\Re }^m\to {\Re }^p$ given by

$r\left(u\right)=y\left(k\right)-\bar{y}\left(k\right),$ (19)

where the model output (18) is reformulated as

$\bar{y}\left(k\right)=C\bar{x}\left(k\right)+Du\left(k\right)$ . (20)

where $D$ is a $p\times m$ coefficient matrix. Formulate Equation (20) as the following input-output equations [23] for $k=0,1,\cdots ,N$ :

$\left[\begin{array}{c}\bar{y}\left(0\right)\\ \bar{y}\left(1\right)\\ \bar{y}\left(2\right)\\ ⋮\\ \bar{y}\left(N\right)\end{array}\right]=\left[\begin{array}{c}C\\ CA\\ C{A}^2\\ ⋮\\ C{A}^N\end{array}\right]{\bar{x}}_0+\left[\begin{array}{ccccc}D& 0& 0& \cdots & 0\\ CB& D& 0& \cdots & 0\\ CAB& CB& D& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ C{A}^{N-1}B& C{A}^{N-2}B& C{A}^{N-3}B& \cdots & D\end{array}\right]\left[\begin{array}{c}u\left(0\right)\\ u\left(1\right)\\ u\left(2\right)\\ ⋮\\ u\left(N-1\right)\end{array}\right].$ (21)

For simplicity, we have

$\bar{y}=E{\bar{x}}_0+Fu$ (22)

where

$E=\left[\begin{array}{c}C\\ CA\\ C{A}^2\\ ⋮\\ C{A}^N\end{array}\right]$ and $F=\left[\begin{array}{ccccc}D& 0& 0& \cdots & 0\\ CB& D& 0& \cdots & 0\\ CAB& CB& D& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ C{A}^{N-1}B& C{A}^{N-2}B& C{A}^{N-3}B& \cdots & D\end{array}\right].$

In addition, define the objective function $g_2:{\Re }^m\to \Re$ , which represents the sum squares of error (SSE), given by

$g_2\left(u\right)=r{\left(u\right)}^{\text{T}}r\left(u\right)$ . (23)

Then, an optimization problem, which is referred to as Problem (O), is stated as follows:

Problem (O):

Find a set of the control sequence $u\left(k\right),k=0,1,\cdots ,N-1$ , such that the objective function $g_2$ is minimized.

It is obviously noticed that for solving Problem (O), Taylor’s theorem [24] [25] is applied to write the objective function $g_2$ as the second-order Taylor expansion about the current $u^{\left(i\right)}$ at iteration i as follows:

$\begin{array}{l}{g}_2\left(u^{\left(i+1\right)}\right)\approx g_2\left(u^{\left(i\right)}\right)+{\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)}^{\text{T}}\nabla g_2\left(u^{\left(i\right)}\right)\\ +\frac{1}{2}{\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)}^{\text{T}}\left({\nabla }^2{g}_2\left(u^{\left(i\right)}\right)\right)\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)\end{array}$ (24)

where the higher-order terms are ignored and the notation $\nabla$ represents the differential operator. The first-order condition in Equation (24) with respect to $u^{\left(i+1\right)}$ is expressed by

$0\approx \nabla g_2\left(u^{\left(i\right)}\right)+\left({\nabla }^2{g}_2\left(u^{\left(i\right)}\right)\right)\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)$ . (25)

Rearrange Equation (25) to yield the normal equation,

$\left({\nabla }^2{g}_2\left(u^{\left(i\right)}\right)\right)\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)=-\nabla g_2\left(u^{\left(i\right)}\right)$ . (26)

Notice that the gradient $\nabla g_2$ is calculated from

$\nabla g_2\left(u^{\left(i\right)}\right)=2\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}r\left(u^{\left(i\right)}\right)$ (27)

and the Hessian of $g_2$ is computed from

${\nabla }^2{g}_2\left(u^{\left(i\right)}\right)=2\left({\nabla }^{2}r{\left(u^{\left(i\right)}\right)}^{\text{T}}r\left(u^{\left(i\right)}\right)+\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}\nabla r\left(u^{\left(i\right)}\right)\right)$ (28)

where $\nabla r\left(u^{\left(i\right)}\right)$ is the Jacobian matrix of $r\left(u^{\left(i\right)}\right)$ , and its entries are denoted by

${\left(\nabla r\left(u^{\left(i\right)}\right)\right)}_{ij}=\frac{\partial r_i}{\partial u_j}\left(u^{\left(i\right)}\right)=F,i=1,2,\cdots ,N-1;j=1,2,\cdots ,N-1.$ (29)

From Equations (27) and (28), Equation (26) can be rewritten as

$\left({\nabla }^{2}r{\left(u^{\left(i\right)}\right)}^{\text{T}}r\left(u^{\left(i\right)}\right)+\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}\nabla r\left(u^{\left(i\right)}\right)\right)\left(u^{\left(i+1\right)}-u^{\left(i\right)}\right)=-\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}r\left(u^{\left(i\right)}\right).$ (30)

By ignoring the second-order derivative term, that is, the first term at the left-hand side of Equation (30), and define

$\Delta u^{\left(i\right)}=u^{\left(i+1\right)}-u^{\left(i\right)}$ , (31)

the normal equation given by Equation (26) is simplified to

$\left(\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}\nabla r\left(u^{\left(i\right)}\right)\right)\Delta u^{\left(i\right)}=-\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}r\left(u^{\left(i\right)}\right)$ . (32)

Then, we obtain the following updating recurrence relation,

$u^{\left(i+1\right)}=u^{\left(i\right)}+\Delta u^{\left(i\right)}$ (33)

with the initial $u^{\left(0\right)}$ given, where

$\Delta u^{\left(i\right)}=-{\left(\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}\nabla r\left(u^{\left(i\right)}\right)\right)}^{-1}\left(\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}r\left(u^{\left(i\right)}\right)\right)$ . (34)

Hence, Equations (33) and (34) are known as the least squares recursive equations, which are based on the Gauss-Newton recursion formula.

From the discussion above, the least-squares updating scheme for the control sequence is summarized below:

Algorithm 2: Least-squares updating scheme

Step 0 Given an initial $u^{\left(0\right)}$ and the tolerance $\epsilon$ . Set $i=0$ .

Step 1 Evaluate the output error $r\left(u^{\left(i\right)}\right)$ and the Jacobian matrix $\nabla r\left(u^{\left(i\right)}\right)$ from Equations (19) and (29), respectively.

Step 2 Solve the normal equation from Equation (32) to obtain $\Delta u^{\left(i\right)}$ .

Step 3 Update the control sequence by using Equation (33). If $u^{\left(i+1\right)}=u^{\left(i\right)}$ , within a given tolerance $\epsilon$ , stop; else set $i=i+1$ and repeat from Step 1 to Step 3.

Remarks:

1) In Step 1, the calculation of the output error $r\left(u^{\left(i\right)}\right)$ and the Jacobian matrix $\nabla r\left(u^{\left(i\right)}\right)$ are done online, however, the Jacobian matrix $\nabla r\left(u^{\left(i\right)}\right)$ could be done offline if it is independent from $u^{\left(i\right)}$ .

2) In Step 2, the inverse of $\nabla r{\left(u^{\left(i\right)}\right)}^{\text{T}}\nabla r\left(u^{\left(i\right)}\right)$ must be exist, and the value of $\Delta u^{\left(i\right)}$ represents the step-size for the control set-point.

3) In Step 3, the initial $u^{\left(0\right)}$ is calculated from (9). The condition $u^{\left(i+1\right)}=u^{\left(i\right)}$ is required to be satisfied for the converged optimal control sequence. The following 2-norm is computed and it is compared with a given tolerance to verify the convergence of $u\left(k\right)$ :

$\Vert u^{\left(i+1\right)}-u^{\left(i\right)}\Vert ={\left({\displaystyle \underset{k=0}{\overset{N-1}{\sum }}\Vert u{\left(k\right)}^{\left(i+1\right)}-u{\left(k\right)}^{\left(i\right)}\Vert }\right)}^{1/2}$ . (35)

4) In order to provide a convergence mechanism for the state sequence, a simple relaxation method is employed:

${\bar{x}}^{\left(i+1\right)}={\bar{x}}^{\left(i\right)}+k_x\left(x-{\bar{x}}^{\left(i\right)}\right)$ (36)

where $k_x\in \left(0,1\right]$ , and $x$ is the state sequence of the real plant.

4. Illustrative Example

In this section, an example of the direct current and alternating current (DC/AC) converter model is illustrated [26] [27] [28] . In this model, the real plant is a nonlinear dynamics structure and the output channel is a single output measurement, where the stationary Gaussian white noises are taken into the model. Since this model is a continuous time model, the simple discretization scheme is required. We define this model as Problem (P), while the simplified model, which is an expectation model, is employed as Problem (M). The cost function for both problems is a quadratic criterion function. The true solution of the original optimal control problem would be obtained by using the approach proposed and the solution procedure is implemented in the MATLAB environment.

Problem (P):

$\underset{u}{\min }{g}_0\left(u\right)=E[\frac{1}{2}x{\left(N\right)}^{\text{T}}S\left(N\right)x\left(N\right)+\frac{1}{2}{\displaystyle \underset{k=0}{\overset{N-1}{\sum }}\left(x{\left(k\right)}^{\text{T}}Qx\left(k\right)+u{\left(k\right)}^{\text{T}}Ru\left(k\right)\right)}]$

subject to

${\stackrel{\dot{}}{x}}_1\left(t\right)=\frac{{\left(x_2\left(t\right)\right)}^2}{x_1\left(t\right)}-5x_1\left(t\right)+5u\left(t\right)+{\omega }_1\left(t\right)$

${\stackrel{\dot{}}{x}}_2\left(t\right)=\frac{{\left(x_2\left(t\right)\right)}^3}{{\left(x_1\left(t\right)\right)}^2}-7x_2\left(t\right)+\left(5\frac{x_2\left(t\right)}{x_1\left(t\right)}+2x_1\left(t\right)\right)u\left(t\right)+{\omega }_2\left(t\right)$

$y\left(t\right)=x_2\left(t\right)+\eta \left(k\right)$

with the initial $x\left(0\right)={\left(0.1,0\right)}^{\text{T}}$ . Here, $\omega \left(k\right)={\left({\omega }_1\left(k\right){\omega }_2\left(k\right)\right)}^{\text{T}}$ and $\eta \left(k\right)$ are Gaussian white noise sequences with their respective covariance given by

$Q_{\omega }={10}^{-2}$ and $R_{\eta }={10}^{-1}$ . The weighting matrices in the cost function are

$S\left(N\right)=I_{2\times 2}$ , $Q=I_{2\times 2}$ and $R=1$ .

Problem (M):

$\underset{u}{\min }{g}_1\left(u\right)=\frac{1}{2}\bar{x}{\left(N\right)}^{\text{T}}S\left(N\right)\bar{x}\left(N\right)+\frac{1}{2}{\displaystyle \underset{k=0}{\overset{N-1}{\sum }}\left(\bar{x}{\left(k\right)}^{\text{T}}Q\bar{x}\left(k\right)+u{\left(k\right)}^{\text{T}}Ru\left(k\right)\right)}$

subject to

$\left(\begin{array}{c}{\bar{x}}_1\left(k+1\right)\\ {\bar{x}}_2\left(k+1\right)\end{array}\right)=\left(\begin{array}{cc}1-5\cdot T& 0\\ 0& 1-7\cdot T\end{array}\right)\left(\begin{array}{c}{\bar{x}}_1\left(k\right)\\ {\bar{x}}_2\left(k\right)\end{array}\right)+\left(\begin{array}{c}5\cdot T\\ 0.2\cdot T\end{array}\right)u\left(k\right)$

$\bar{y}\left(k\right)={\bar{x}}_2\left(k\right)$

for $k=0,1,\cdots ,80$ , with the sampling time $T=0.01$ second.

The simulation result is shown in Table 1. The implementation of the algorithm proposed takes four iteration numbers to converge. The initial cost of 0.0429 unit is the value of cost function belong to Problem (M) before the iteration begins. After the convergence is achieved, the final cost of 116.9461 units is obtained. The original cost of 1.0881 × 10³ is the value of cost function of the original optimal control problem. It is noticed that the cost reduction is 89.28 percent. The value of SSE, which is 8.046097 × 10⁻¹², shows the smallest differences between the real output and the model output.

The trajectories of final control and real control are shown in Figure 1 and Figure 2, respectively. The final control takes into the consideration of noise sequence during the calculation, while the real control is free from the noise disturbance. Figure 3 shows the trajectories of final output and real output. It can be seen that the final output tracks closely to the real output. Figure 4 shows the trajectories of final state and real state, where they are statistically identical. Figure 5 and Figure 6 show, respectively, the output errors after and before iterations. The accuracy of the model output is increased due on the decrease of the output error.

5. Concluding Remarks

A computational algorithm, which is equipped with the efficient matching scheme, was discussed in this paper. The linear model-based optimal control problem, which is simplified from the discrete time nonlinear stochastic optimal control problem, was solved iteratively. During the calculation procedure, the model used was reformulated into the input-output equations. Then, the least squares optimization problem was introduced. By satisfying the first order necessary condition, the normal equation was solved such that the least squares

recursion formula could be established. In this way, the control policy could be updated iteratively. As a result of this, the sum squares of the output errors was minimized, which indicates the equivalent of the model output to the real output. For illustration, an example of the direct current and alternating current converter model was studied. The result obtained showed the accuracy of the algorithm proposed. In conclusion, the efficiency of the algorithm proposed was highly presented.

Acknowledgements

The authors would like to thanks the Universiti Tun Hussein Onn Malaysia (UTHM) for financial supporting to this study under Incentive Grant Scheme for Publication (IGSP) VOT. U417. The second author was supported by the NSF (11501053) of China and the fund (15C0026) of the Education Department of Hunan Province.

Cite this paper

Kek, S.L., Li, J. and Teo, K.L. (2017) Least Squares Solution for Discrete Time Nonlinear Stochastic Optimal Control Problem with Model-Reality Differences. Applied Mathematics, 8, 1-14. http://dx.doi.org/10.4236/am.2017.81001

References