Conjugate Gradient Approach for Discrete Time Optimal Control Problems with Model-Reality Differences

Abstract

In this chapter, an efficient computation approach is proposed for solving a general class of discrete-time optimal control problems. In our approach, a simplified optimal control model, which is adding the adjusted parameters into the model used, is solved iteratively. In this way, the differences between the real plant and the model used are calculated, in turn, to update the optimal solution of the model used. During the computation procedure, the equivalent optimization problem is formulated, where the conjugate gradient algorithm is applied in solving the optimization problem. On this basis, the optimal solution of the modified model-based optimal control problem is obtained repeatedly. Once the convergence is achieved, the iterative solution approximates to the correct optimal solution of the original optimal control problem, in spite of model-reality differences. For illustration, both linear and nonlinear examples are demonstrated to show the performance of the approach proposed. In conclusion, the efficiency of the approach proposed is highly presented.

1. Introduction

Optimal control problems are existing in engineering and natural sciences for so long, and the applications of the optimal control have been well defined in the literature [1, 2, 3, 4]. With the rapid evolution of computer technology, the development of the optimal control techniques is reached a mature level, from classical control to modern control, from proportional-integral-derivative (PID) control to feedback control, and from adaptive control to intelligent control [5, 6, 7, 8]. The studies in the optimal control field are still progressing, and attract the interest of, not only engineers and applied mathematicians but also biologists and financialists, to investigate and contribute to the optimal control theory.

In particular, the optimal control algorithm, which integrates system optimization and parameter estimation, gives a new insight into the control community. This algorithm is known as the integrated system optimization and parameter estimation (ISOPE), and its dynamic version is called the dynamic ISOPE (DISOPE). Both of these algorithms were introduced by Robert [9, 10, 11], and Robert and Becerra [12, 13, 14], respectively. The basic idea of DISOPE is applying the model-based optimal control, which has different structures and parameters compared to the original optimal control problem, to obtain the correct optimal solution of the original optimal control problem, in spite of model-reality differences. Recently, this algorithm has been extended to cover both deterministic and stochastic versions, and it is known as an integrated optimal control and parameter estimation (IOCPE) algorithm [15, 16]. On the other hand, the application of the optimization techniques, particularly, using the conjugate gradient method for solving the optimal control problem [17, 18, 19] has also been studied, where the open-loop control strategy is concerned [3, 8].

In this chapter, the conjugate gradient approach [17, 19] is employed to solve the linear model-based optimal control problem for obtaining the optimal solution of the original optimal control problem. In our approach, the simplified model, which is adding the adjusted parameters, is formulated initially. Then, an expanded optimal control problem, which combines the system dynamic and the cost function from the original optimal control problem into the simplified model, is introduced. By defining the Hamiltonian function and the augmented cost function, the corresponding necessary conditions for optimality are derived. Among these necessary conditions, a set of necessary conditions is for the modified model-based optimal control problem, a set of necessary conditions defines the parameter estimation problem, and a set of necessary conditions calculates the multipliers [15].

By virtue of the modified model-based optimal control problem, an equivalence optimization problem is defined, and the related gradient function is determined. With an initial control sequence, the initial gradient and the initial search direction are computed. Then, the control sequences are updated through the line search technique, where the gradient and search direction would satisfy the conjugacy condition [17, 18]. During the iteration, the state and the costate are updated by the control sequence obtained from the conjugate gradient approach. When the convergence is achieved within a tolerance given, the iterative solution approximates to the correct optimal solution of the original optimal control problem, in spite of model-reality differences. For illustration, examples of linear and nonlinear cases, which are damped harmonic oscillator [7] and continuous stirred-tank chemical tank [8], are studied.

The chapter is organized as follows. In Section 2, the problem statement is described in detail, where the original optimal control problem and the simplified model are discussed. In Section 3, the methodology used is further explained. The necessary conditions for optimality are derived, and the use of the conjugate gradient method is delivered in solving the equivalence optimization problem. In Section 4, examples of a damped harmonic oscillator and a continuous stirred-tank chemical reactor are studied. The results show the efficiency of the algorithm proposed. Finally, concluding remarks are made.

2. Problem statement

Consider a general class of the discrete-time nonlinear optimal control problem, given by

$$\begin{array}{c}\underset{u\left(k\right)}{min}{J}_0\left(u\right)=\phi \left(x\left(N\right),N\right)+\sum _{k=0}^{N-1}L\left(x\left(k\right),u\left(k\right),k\right)\\ \text{subject to}\hfill \\ x\left(k+1\right)=f\left(x\left(k\right),u\left(k\right),k\right),x\left(0\right)=x_0\end{array}$$ E1

where $u\left(k\right)\in {\Re }^m,$ $k=0,1,\cdots ,N-1,$ and $x\left(k\right)\in {\Re }^n,$ $k=0,1,\cdots ,N,$ are the control sequences and the state sequences, respectively, while $f:{\Re }^n\times {\Re }^m\times \Re \to {\Re }^n$ represents the real plant, $L:{\Re }^n\times {\Re }^m\times \Re \to \Re$ is the cost under summation, and $\phi :{\Re }^n\times \Re \to \Re$ is the terminal cost. Here, $J_0$ is the scalar cost function, and $x_0$ is the known initial state vector. It is assumed that all functions in (1) are continuously differentiable with respect to their respective arguments.

This problem, which is referred to as Problem (P), is regarded as the real optimal control problem. Due on the complex and nonlinear structure, solving Problem (P) actually requires the efficient computation techniques. For this reason, the simplified model of Problem (P) is identified to be solved such that the true optimal solution of Problem (P) could be approximated. Hence, this simplified model-based optimal control problem is defined as follows:

$$\begin{array}{l}\underset{u\left(k\right)}{min}{J}_1\left(u\right)=\frac{1}{2}x{\left(N\right)}^{\mathrm{T}}S\left(N\right)x\left(N\right)+\gamma \left(N\right)+\sum _{k=0}^{N-1}\frac{1}{2}\left(x{\left(k\right)}^{\mathrm{T}}\mathit{Qx}\left(k\right)+u{\left(k\right)}^{\mathrm{T}}\mathit{Ru}\left(k\right)\right)+\gamma \left(k\right)\\ \text{subject to}\\ x\left(k+1\right)=\mathit{Ax}\left(k\right)+\mathit{Bu}\left(k\right)+\alpha \left(k\right),x\left(0\right)=x_0\end{array}$$ E2

where $\gamma \left(k\right),k=0,1,\cdots ,N,$ and $\alpha \left(k\right),k=0,1,\cdots ,N-1,$ are introduced as the adjusted parameters, whereas $A$ is an $n\times n$ transition matrix, and $B$ is an $n\times m$ control coefficient matrix. Besides, $S\left(N\right)$ and $Q$ are $n\times n$ positive semi-definite matrices, and $R$ is a $m\times m$ positive definite matrix. Here, $J_1$ is the scalar cost function.

Let this problem is referred to as Problem (M). It can be seen that, because of the different structures and parameters, only solving Problem (M) would not obtain the optimal solution of Problem (P) for not taking the adjusted parameters into account. Notice, adding the adjusted parameters into Problem (M) could let us calculate the differences between the real plant and the model used. On this basis, Problem (M) would be solved iteratively to give the correct optimal solution of Problem (P), in spite of model-reality differences.

3. System optimization with parameter estimation

Now, an expanded optimal control problem, which combines the real plant and the cost function in Problem (P) into Problem (M) and is referred to as Problem (E), is introduced by

$$\begin{array}{l}\underset{u\left(k\right)}{min}{J}_2\left(u\right)=\frac{1}{2}x{\left(N\right)}^{\mathrm{T}}S\left(N\right)x\left(N\right)+\gamma \left(N\right)+\sum _{k=0}^{N-1}\frac{1}{2}\left(x{\left(k\right)}^{\mathrm{T}}\mathit{Qx}\left(k\right)+u{\left(k\right)}^{\mathrm{T}}\mathit{Ru}\left(k\right)\right)+\gamma \left(k\right)\\ +\frac{1}{2}{r}_1{‖u\left(k\right)-v\left(k\right)‖}^2+\frac{1}{2}{r}_2{‖x\left(k\right)-z\left(k\right)‖}^2\\ \text{subject to}\\ x\left(k+1\right)=\mathit{Ax}\left(k\right)+\mathit{Bu}\left(k\right)+\alpha \left(k\right),x\left(0\right)=x_0\\ \frac{1}{2}z{\left(N\right)}^{\mathrm{T}}S\left(N\right)z\left(N\right)+\gamma \left(N\right)=\phi \left(z\left(N\right),N\right)\\ \frac{1}{2}\left(z{\left(k\right)}^{\mathrm{T}}\mathit{Qz}\left(k\right)+v{\left(k\right)}^{\mathrm{T}}\mathit{Rv}\left(k\right)\right)+\gamma \left(k\right)=L\left(z\left(k\right),v\left(k\right),k\right)\\ \mathit{Az}\left(k\right)+\mathit{Bv}\left(k\right)+\alpha \left(k\right)=f\left(z\left(k\right),v\left(k\right),k\right)\\ v\left(k\right)=u\left(k\right)\\ z\left(k\right)=x\left(k\right)\end{array}$$ E3

where $v\left(k\right)\in {\Re }^m,$ $k=0,1,\cdots ,N-1,$ and $z\left(k\right)\in {\Re }^n,$ $k=0,1,\cdots ,N,$ are introduced to separate the sequences of control and state in the optimization problem from the respective signals in the parameter estimation problem, and $\Vert \cdot \Vert$ denotes the usual Euclidean norm. The terms $\frac{1}{2}{r}_1{‖u\left(k\right)-v\left(k\right)‖}^2$ and $\frac{1}{2}{r}_2{‖x\left(k\right)-z\left(k\right)‖}^2$ with $r_1,r_2\in \Re$ are introduced to improve the convexity and to facilitate the convergence of the resulting iterative algorithm. Here, we clarify that the algorithm is designed such that the constraints $v\left(k\right)=u\left(k\right)$ and $z\left(k\right)=x\left(k\right)$ are satisfied upon termination of the iterations, assuming that convergence is achieved. Moreover, the state constraint $z\left(k\right)$ and the control constraint $v\left(k\right)$ are used for the computation of the parameter estimation and matching scheme, while the corresponding state constraint $x\left(k\right)$ and control constraint $u\left(k\right)$ are reserved for optimizing the model-based optimal control problem. Therefore, system optimization and parameter estimation are declared and mutually integrated.

3.1 Necessary conditions for optimality

Define the Hamiltonian function for Problem (E), given by:

$$\begin{array}{ll}{H}_2(k)=& \frac{1}{2}(x{(k)}^{\text{T}}Qx(k)+u{(k)}^{\text{T}}Ru(k))+\gamma (k)\\ & +\frac{1}{2}{r}_1\left|\right|u(k)-v(k)|{|}^2+\frac{1}{2}{r}_2||x(k)-z(k)|{|}^2\\ & +p{(k+1)}^{\text{T}}(Ax(k)+Bu(k)+\alpha (k))\\ & -\lambda {(k)}^{\text{T}}u(k)-\beta {(k)}^{\text{T}}x(k)\end{array}$$ E4

where $\lambda \left(k\right)\in {\Re }^m,$ $k=0,1,\cdots ,N-1,\beta \left(k\right)\in {\Re }^n,$ $k=0,1,\cdots ,N,$ and $p\left(k\right)\in {\Re }^n,k=0,1,\cdots ,N,$ are modifiers. Using this Hamiltonian function in (4), write the cost function in (3) to be the augmented cost function, that is,

$$\begin{array}{ll}{J^{\prime }}_2(u)=& \frac{1}{2}x{(N)}^{\text{T}}S(N)x(N)+\gamma (N)+p{(0)}^{\text{T}}x(0)-p{(N)}^{\text{T}}x(N)\\ & +\xi (N)(\phi (z(N),N)-\frac{1}{2}z{(N)}^{\text{T}}S(N)z(N)-\gamma (N))\\ & +{\Gamma }^{\text{T}}(x(N)-z(N))+{\displaystyle \sum _{k=0}^{N-1}{H}_2(k)}-p{(k)}^{\text{T}}x(k)+\lambda {(k)}^{\text{T}}v(k)+\beta {(k)}^{\text{T}}z(k)\\ & +\xi (k)(L(z(k),v(k),k)-\frac{1}{2}(z{(k)}^{\text{T}}Qz(k)+v{(k)}^{\text{T}}Rv(k))-\gamma (k))\\ & +\mu {(k)}^{\text{T}}(f(z(k),v(k),k)-Az(k)-Bv(k)-\alpha (k))\end{array}$$ E5

where $p\left(k\right),$ $\xi \left(k\right),$ $\lambda \left(k\right),$ $\beta \left(k\right),$ $\mu \left(k\right)$ and $\mathrm{\Gamma }$ are the appropriate multipliers to be determined later.

Applying the calculus of variation [7, 9, 11, 13, 15] to the augmented cost function in (5), the following necessary conditions for optimality are obtained:

(a) Stationary condition:

$$\mathit{Ru}\left(k\right)+B^{\mathrm{T}}p\left(k+1\right)-\lambda \left(k\right)+r_1\left(u\left(k\right)-v\left(k\right)\right)=0$$ E6

(b) Co-state equation:

$$p\left(k\right)=\mathit{Qx}\left(k\right)+A^{\mathrm{T}}p\left(k+1\right)-\beta \left(k\right)+r_2\left(x\left(k\right)-z\left(k\right)\right)$$ E7

(c) State equation:

$$x\left(k+1\right)=\mathit{Ax}\left(k\right)+\mathit{Bu}\left(k\right)+\alpha \left(k\right)$$ E8

(d) Boundary conditions:

$$p\left(N\right)=S\left(N\right)x\left(N\right)+\mathrm{\Gamma }\text{and}x\left(0\right)=x_0$$ E9

(e) Adjusted parameter equations:

$$\phi \left(z\left(N\right),N\right)=\frac{1}{2}z{\left(N\right)}^{\mathrm{T}}S\left(N\right)z\left(N\right)+\gamma \left(N\right)$$ E10

$$L\left(z\left(k\right),v\left(k\right),k\right)=\frac{1}{2}\left(z{\left(k\right)}^{\mathrm{T}}\mathit{Qz}\left(k\right)+v{\left(k\right)}^{\mathrm{T}}\mathit{Rv}\left(k\right)\right)+\gamma \left(k\right)$$ E11

$$f\left(z\left(k\right),v\left(k\right),k\right)=\mathit{Az}\left(k\right)+\mathit{Bv}\left(k\right)+\alpha \left(k\right)$$ E12

(f) Modifier equations:

$$\mathrm{\Gamma }={\nabla }_{z\left(N\right)}\phi -S\left(N\right)z\left(N\right)$$E13

$$\lambda \left(k\right)=-\left({\nabla }_{v\left(k\right)}L-\mathit{Rv}\left(k\right)\right)-{\left(\frac{\partial f}{\partial v\left(k\right)}-B\right)}^{\mathrm{T}}\widehat{p}\left(k+1\right)$$ E14

$$\beta \left(k\right)=-\left({\nabla }_{z\left(k\right)}L-\mathit{Qz}\left(k\right)\right)-{\left(\frac{\partial f}{\partial z\left(k\right)}-A\right)}^{\mathrm{T}}\widehat{p}\left(k+1\right)$$ E15

with $\xi \left(k\right)=1$ and $\mu \left(k\right)=\widehat{p}\left(k+1\right)$ .

(g) Separable variables:

$$v\left(k\right)=u\left(k\right),z\left(k\right)=x\left(k\right),\widehat{p}\left(k\right)=p\left(k\right).$$ E16

Notice that for the optimality necessary conditions obtained above, they are divided into three sets of necessary conditions. The first set of necessary conditions in (6)–(9) is the necessary conditions for the system optimization problem. The second set of necessary conditions in (10)–(12) defines the parameter estimation problem. The third set of necessary conditions in (13)–(15) provides the computation of multipliers. In fact, the necessary conditions, which are defined in (6)–(9), are the optimality for the modified model-based optimal control problem, and the adjusted parameters, which are calculated from the necessary conditions in (10)–(12), measure the differences between the real plant and the model used.

3.2 Modified model-based optimal control problem

As a consequence, the modified model-based optimal control problem, which is referred to as Problem (MM), is defined by

$$\begin{array}{l}\underset{u\left(k\right)}{min}{J}_3\left(u\right)=\frac{1}{2}x{\left(N\right)}^{\mathrm{T}}S\left(N\right)x\left(N\right)+{\mathrm{\Gamma }}^{\mathrm{T}}x\left(N\right)+\gamma \left(N\right)+\sum _{k=0}^{N-1}\frac{1}{2}\left(x{\left(k\right)}^{\mathrm{T}}\mathit{Qx}\left(k\right)+u{\left(k\right)}^{\mathrm{T}}\mathit{Ru}\left(k\right)\right)\\ +\gamma \left(k\right)+\frac{1}{2}{r}_1{‖u\left(k\right)-v\left(k\right)‖}^2+\frac{1}{2}{r}_2{‖x\left(k\right)-z\left(k\right)‖}^2-\lambda {\left(k\right)}^{\mathrm{T}}u\left(k\right)-\beta {\left(k\right)}^{\mathrm{T}}x\left(k\right)\\ \text{subject to}\\ x\left(k+1\right)=\mathit{Ax}\left(k\right)+\mathit{Bu}\left(k\right)+\alpha \left(k\right),x\left(0\right)=x_0\end{array}$$ E17

with the specified $\alpha \left(k\right),$ $\gamma \left(k\right),$ $\lambda \left(k\right),$ $\beta \left(k\right),$ $\mathrm{\Gamma },$ $v\left(k\right)$ and $z\left(k\right),$ where the boundary conditions are given by $x_0$ and $p\left(N\right)$ with the specified multiplier $\mathrm{\Gamma }$ .

It is obvious that Problem (MM), which is derived from Problem (E), is a modification of optimal control problem and is also known as a modified linear quadratic regular problem. Importantly, the set of the necessary conditions in (6)–(9) for Problem (E) is the necessary conditions that are satisfied by Problem (MM). In addition, due to the quadratic criterion feature of the objective function, the conjugate gradient method [17, 18], which is one of the numerical optimization techniques, could be applied to solve Problem (MM).

3.3 Conjugate gradient algorithm

For simplicity [19], establish Problem (MM) as a nonlinear optimization problem with the initial control given by $u^{\left(0\right)}=u{\left(k\right)}^0$ as follows:

$$\underset{u\left(k\right)}{min}{J}_3\left(u\right)\text{subject to}u=u\left(k\right)\in {\Re }^m\text{for}k=0,1,\cdots ,N-1$$E18

Let this problem as Problem (Q). Moreover, the Hamiltonian function defined in (4) is taken into consideration as an equivalent objective function. Hence, this Hamiltonian function allows the evaluation of the gradient function, which is the stationary condition in (6), and by using the iterative solution $u^{\left(i\right)}=u{\left(k\right)}^i$ to satisfy the state Eq. (8), which is solved forward in time, and the co-state Eq. (7), which is solved backward in time.

Define the gradient function $g:{\Re }^m\to {\Re }^m$ as

$$g\left(u^i\right)={\nabla }_u{J}_3\left(u^i\right)$$ E19

which is represented by the stationary condition in (6). For arbitrary initial control $u^{\left(0\right)}\in {\Re }^m$ , the initial gradient and the initial search direction are calculated from

$$g^{\left(0\right)}=g\left(u^{\left(0\right)}\right)$$ E20

$$d^{\left(0\right)}=-g^{\left(0\right)}.$$E21

The following line search equation is applied to update the control sequence:

$$u^{\left(i+1\right)}=u^{\left(i\right)}+a_i\cdot d^{\left(i\right)}$$ E22

where $a_i\in \Re$ is the step size, and its value can be determined from

$$a_i=\underset{a\ge 0}{argmin}{J}_3\left(u^{\left(i\right)}+a\cdot d^{\left(i\right)}\right).$$E23

After that, the gradient and the search direction are updated by

$$g^{\left(i+1\right)}=g\left(u^{\left(i+1\right)}\right)$$ E24

$$d^{\left(i+1\right)}=-g^{\left(i+1\right)}+b_i\cdot d^{\left(i\right)}$$ E25

with

$$b_i=\frac{g^{\left(i+1\right)\mathrm{T}}\cdot g^{\left(i+1\right)}}{g^{\left(i\right)\mathrm{T}}\cdot g^{\left(i\right)}}$$ E26

for i = 0, 1, 2, … represents the iteration numbers.

From the discussion above, we present the result as a proposition given as follows:

Proposition 1. Consider Problem (Q). The control sequence $u^{\left(i\right)}$ , which is defined in (22) and is represented by

$$u^{\left(i\right)}=\left[{\left(u\left(0\right)\right)}^{\mathrm{T}},{\left(u\left(1\right)\right)}^{\mathrm{T}},\cdots ,{\left(u\left(N-1\right)\right)}^{\mathrm{T}}\right],$$

is generated through a set of the search direction vector $d^{\left(i\right)}$ whose components are linearly independent. Also, the direction $d^{\left(i\right)}$ is conjugacy.

The conjugate gradient algorithm is summarized below:

Conjugate gradient algorithm

Data: Choose the arbitrary initial control $u^{\left(0\right)}$ and the tolerance $\epsilon$ .

Step 0: Compute the initial gradient $g^{\left(0\right)}$ from (20) and the initial search direction $d^{\left(0\right)}$ from (21), respectively. Set $i=0.$

Step 1: Solve the state Eq. (8) forward in time from $k=0$ to $k=N$ with the initial condition (9) to obtain $x{\left(k\right)}^i,k=0,1,\cdots ,N.$

Step 2: Solve the costate Eq. (7) backward in time from $k=N$ to $k=0$ with the boundary condition (9), where $p{\left(k\right)}^i$ is the solution obtained.

Step 3: Calculate the value of the cost functional $J_3\left(u^{\left(i\right)}\right)$ from (17).

Step 4: Solve (23) to obtain the step size $a_i$ .

Step 5: Calculate the control $u^{\left(i+1\right)}$ from (22).

Step 6: Evaluate the gradient $g^{\left(i+1\right)}$ and the search direction $d^{\left(i+1\right)},$ respectively, from (24) and (25) with computing $b_i$ from (26). If the gradient $g^{\left(i+1\right)}=g^{\left(i\right)},$ within a given tolerance, stop, else set $i=i+1,$ go to Step 1.

Remark 1:

1. Step 0 is the preliminary step for setting the initial search direction based on the gradient direction in using the conjugate gradient algorithm.

2. Steps 1, 2, and 3 are performed to solve the system optimization by using the corresponding control sequence $u^{\left(i\right)}$ .

3. Steps 4, 5, and 6 are the computation steps in implementing the conjugate direction.

3.4 Iterative calculation procedure

Accordingly, Problem (Q) is solved by using the conjugate gradient algorithm. Indeed, the solution procedure for system optimization with parameter estimation could be described by joining the conjugate gradient algorithm with the parameters estimated. A summary of the calculation procedure including the principle of model-reality differences is listed as follows:

Iterative algorithm based on model-reality differences

Data: $A,$ $B,$ $Q,$ $R,$ $S\left(N\right),$ $N,$ $x_0,$ $r_1,$ $r_2,$ $k_v,$ $k_z,$ $k_p,$ $f,$ $L.$ Note that $A$ and $B$ could be determined based on the linearization of $f$ at $x_0$ or from the linear terms of $f$ .

Step 0: Compute a nominal solution. Assume that $\alpha \left(k\right)=0,$ $k=0,1,\cdots ,N-1,$ and $r_1=r_2=0.$ Solve Problem (M) defined by (2) to obtain $u{\left(k\right)}^0,$ $k=0,1,\cdots ,N-1,$ and $x{\left(k\right)}^0,$ $p{\left(k\right)}^0,$ $k=0,1,\cdots ,N.$ Then, with $\alpha \left(k\right)=0,k=0,1,\cdots ,N-1,$ and using $r_1,$ $r_2$ from the data. Set $i=0,$ $v{\left(k\right)}^0=u{\left(k\right)}^0,$ $z{\left(k\right)}^0=x{\left(k\right)}^0$ and $\widehat{p}{\left(k\right)}^0=p{\left(k\right)}^0.$

Step 1: Compute the parameters $\gamma {\left(k\right)}^i,$ $k=0,1,\cdots ,N,$ and $\alpha {\left(k\right)}^i,$ $k=0,1,\cdots ,N-1,$ from (10)–(12). This is called the parameter estimation step.

Step 2: Compute the modifiers ${\mathrm{\Gamma }}^i,$ $\lambda {\left(k\right)}^i$ and $\beta {\left(k\right)}^i,$ $k=0,1,\cdots ,N-1,$ from (13)–(15). Notice that this step requires taking the derivatives of $f$ and $L$ with respect to $v{\left(k\right)}^i$ and $z{\left(k\right)}^i.$

Step 3: With $\gamma {\left(k\right)}^i,$ $\alpha {\left(k\right)}^i,$ ${\mathrm{\Gamma }}^i,\lambda {\left(k\right)}^i,$ $\beta {\left(k\right)}^i,$ $v{\left(k\right)}^i$ , and $z{\left(k\right)}^i,$ solve Problem (Q) using the conjugate gradient algorithm. This is called the system optimization step.

Step 4: Test the convergence and update the optimal solution of Problem (P). In order to provide a mechanism for regulating convergence, a simple relaxation method is employed:

$$v{\left(k\right)}^{i+1}=v{\left(k\right)}^i+k_v\left(u{\left(k\right)}^i-v{\left(k\right)}^i\right)$$ E27

$$z{\left(k\right)}^{i+1}=z{\left(k\right)}^i+k_z\left(x{\left(k\right)}^i-z{\left(k\right)}^i\right)$$ E28

$$\widehat{p}{\left(k\right)}^{i+1}=\widehat{p}{\left(k\right)}^i+k_p\left(p{\left(k\right)}^i-\widehat{p}{\left(k\right)}^i\right)$$ E29

where $k_v,$ $k_z,$ $k_p\in \left(0,1\right]$ are scalar gains. If $v{\left(k\right)}^{i+1}=v{\left(k\right)}^i,$ $k=0,1,\cdots ,N-1,$ and $z{\left(k\right)}^{i+1}=z{\left(k\right)}^i,$ $k=0,1,\cdots ,N,$ within a given tolerance, stop; else set $i=i+1,$ and repeat the procedure starting from Step 1.

Remark 2:

1. In Step 0, the nominal solution could be obtained by using the standard procedure of the linear quadratic regulator approach, where the feedback gain and the Riccati equation are calculated offline.

2. In Step 3, applying the conjugate gradient algorithm to obtain the new control sequence will give a good effect if the conjugacy of the search direction is satisfied.

3. In Step 4, the simple relaxation method in (27)–(29) is used, so that the matching scheme for the parameters and the optimal solution can be established.

4. Illustrative examples

In this section, two examples are studied. The first example is for optimizing and controlling a damped harmonic oscillator [7], and the second example is related to optimal control of a continuous stirred-tank chemical reactor [8]. The mathematical models of these examples are discussed, and their optimal solution is obtained by using the algorithm discussed in Section 3. Here, the algorithm is implemented in the Octave 5.1.0 environment.

4.1 Example 1: a damped harmonic oscillator

Consider a damped harmonic oscillator [7] given by

$$\dot{x}=\left(\begin{array}{cc}0& 1\\ -{\omega }^2& -2\mathit{\delta \omega }\end{array}\right)x+\left(\begin{array}{c}0\\ 1\end{array}\right)u$$ E30

with the natural frequency $\omega$ = 0.8, the damping ratio $\delta$ = 0.1, and the initial state $x_0={\left(1010\right)}^{\mathrm{T}}.$ Define the state $x={\left(\begin{array}{cc}{x}_1& x_2\end{array}\right)}^{\mathrm{T}},$ where $x_1$ is the displacement and $x_2$ is the velocity. For the purpose of controlling this oscillator, the following objective function

$$J\left(0\right)=\frac{1}{2}{\int }_{0.0}^{9.4}\left({\left(x_1\left(t\right)\right)}^2+{\left(x_2\left(t\right)\right)}^2+{\left(u\left(t\right)\right)}^2\right)\mathit{dt}$$ E31

is minimized. This problem is a continuous-time linear optimal control problem, and the equivalence discrete time optimal control problem, which is regarded as Problem (P), is given by:

$$\begin{array}{l}\underset{u}{min}J\left(u\right)=\sum _{k=0}^{10}\frac{1}{2}\mathrm{\Delta }t\left(x_1{\left(k\right)}^2+x_2{\left(k\right)}^2+u{\left(k\right)}^2\right)\\ \text{subject to}\\ x\left(k+1\right)=\left(\begin{array}{cc}1.00& 0.94\\ -0.60& 0.85\end{array}\right)x\left(k\right)+\left(\begin{array}{c}0.00\\ 0.94\end{array}\right)u\left(k\right)\end{array}$$ E32

with the initial state $x_0={\left(\begin{array}{cc}10& 10\end{array}\right)}^{\mathrm{T}}.$ and the sampling time $\mathrm{\Delta }t=0.94$ s is taken for the discretization transform.

Consider the model-based optimal control problem, which is regarded as Problem (M), given by:

$$\begin{array}{l}\underset{u}{min}J\left(u\right)=\sum _{k=0}^{10}\left(\frac{1}{2}\left(x_1{\left(k\right)}^2+x_2{\left(k\right)}^2+u{\left(k\right)}^2\right)+\gamma \left(k\right)\right)\mathrm{\Delta }t\\ \text{subject to}\\ x\left(k+1\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)x\left(k\right)+\left(\begin{array}{c}1\\ 0\end{array}\right)u\left(k\right)+\alpha \left(k\right)\end{array}$$ E33

with the initial state $x_0={\left(\begin{array}{cc}10& 10\end{array}\right)}^{\mathrm{T}},$ and the adjusted parameters $\gamma \left(k\right)$ , $k=0,1,\cdots ,N,$ and $\alpha \left(k\right)$ , $k=0,1,\cdots ,N-1,$ are supplied to the model used.

By using the algorithm proposed, the simulation result is shown in Table 1. Notice that the minimum cost for Problem (M) is 546.05 units without adding the adjusted parameters. Once the adjusted parameters are taken into consideration, the iterative solution approximates to the true optimal solution of the original optimal control problem, in spite of model-reality differences. It is highlighted that there is a 99% of the cost reduction to obtain the final cost of 128.50 units.

Number of iteration	Initial cost	Final cost	Elapsed time (s)
20	17053.11	128.50	1.38021

Table 1.

Simulation result, Example 1.

Figures 1 and 2 show the trajectories of control and state, respectively. With this control effort, the state reaches at the steady state after 4 units of time, which presents the oscillator stopped from moving. Figure 3 shows the changes of the costate at the first 2 units of time. The optimal solution obtained is verified by satisfying the stationary condition as shown in Figure 4. Figures 5 and 6 show the adjusted parameters after the convergence is achieved, where the model-reality differences are measured during the iterative procedure.

Figure 1.

Final control u(k), Example 1.

Figure 2.

Final state x(k), Example 1.

Figure 3.

Final costate p(k), Example 1.

Figure 4.

Stationary H_u(k), Example 1.

Figure 5.

Adjusted parameter α(k), Example 1.

Figure 6.

Adjusted parameter γ(k), Example 1.

Therefore, this damped harmonic oscillator is controlled, and the cost function is minimized as desired.

4.2 Example 2: a continuous stirred-tank chemical reactor

Consider a continuous stirred-tank chemical reactor, which consists of two state equations [8]. The flow of a coolant through a coil inserted in the reactor is to control the first order, irreversible exothermic reaction taking place in the reactor. Assume that $x_1\left(t\right)$ is the deviation from the steady-state temperature, $x_2\left(t\right)$ is the deviation from the steady-state concentration, and $u\left(t\right)$ is the normalized control variable that represents the effect of coolant flow on the chemical reaction. The corresponding state equations are given by

$${\dot{x}}_1\left(t\right)=-2\left(x_1\left(t\right)+0.25\right)+\left(x_2\left(t\right)+0.5\right)exp\left(\frac{25x_1\left(t\right)}{x_1\left(t\right)+2}\right)-\left(x_1\left(t\right)+0.25\right)u\left(t\right)$$E34

$${\dot{x}}_2\left(t\right)=0.5-x_2\left(t\right)-\left(x_2\left(t\right)+0.5\right)exp\left(\frac{25x_1\left(t\right)}{x_1\left(t\right)+2}\right)$$E35

with the initial state $x_0={\left(\begin{array}{cc}0.05& 0.00\end{array}\right)}^{\mathrm{T}}$ . The cost function to be minimized is given by

$$J\left(0\right)={\int }_{0.0}^{0.8}\left({\left(x_1\left(t\right)\right)}^2+{\left(x_2\left(t\right)\right)}^2+0.1{\left(u\left(t\right)\right)}^2\right)\mathit{dt}.$$E36

Here, the desired objective is to maintain the temperature and the concentration close to their respective steady-state values without expending large amounts of the control effort.

This problem is a continuous time nonlinear optimal control problem. For doing the discretization transform, the sampling time $\mathrm{\Delta }t=0.0057$ s is used to formulate the equivalence discrete-time optimal control problem, which is referred to as Problem (P), given by:

$$\begin{array}{l}\underset{u}{min}J\left(u\right)=\sum _{k=0}^{80}\frac{1}{2}\left(2x_1{\left(k\right)}^2+2x_2{\left(k\right)}^2+0.2u{\left(k\right)}^2\right)\mathrm{\Delta }t\\ \text{subject to}\\ x_1\left(k+1\right)=x_1\left(k\right)-2\left(x_1\left(k\right)+0.25\right)\mathrm{\Delta }t+\left(x_2\left(k\right)+0.5\right)\mathrm{\Delta }texp\left(\frac{25x_1\left(k\right)}{x_1\left(k\right)+2}\right)\\ -\left(x_1\left(k\right)+0.25\right)u\left(k\right)\mathrm{\Delta }t)\\ x_2\left(k+1\right)=x_2\left(k\right)+\left(0.5-x_2\left(k\right)\right)\mathrm{\Delta }t-\left(x_2\left(k\right)+0.5\right)\mathrm{\Delta }texp\left(\frac{25x_1\left(k\right)}{x_1\left(k\right)+2}\right)\end{array}$$ E37

with the initial state $x_0={\left(\begin{array}{cc}0.05& 0.00\end{array}\right)}^{\mathrm{T}}$ .

By applying the algorithm proposed to obtain the optimal solution for Problem (P), the following model, which is referred to as Problem (M), is introduced,

$$\begin{array}{c}\underset{u}{min}J\left(u\right)=\sum _{k=0}^{80}\frac{1}{2}\left(2x_1{\left(k\right)}^2+2x_2{\left(k\right)}^2+0.2u{\left(k\right)}^2\right)\mathrm{\Delta }t\\ \text{subject to}\\ x\left(k+1\right)=\left(\begin{array}{cc}1.048& 0.010\\ -0.062& 0.984\end{array}\right)x\left(k\right)+\left(\begin{array}{c}-0.002\\ 0.000\end{array}\right)u\left(k\right)+\alpha \left(k\right)\end{array}$$ E38

with the initial state $x_0={\left(\begin{array}{cc}0.05& 0.00\end{array}\right)}^{\mathrm{T}},$ and the adjusted parameters $\gamma \left(k\right)$ , $k=0,1,\cdots ,N,$ and $\alpha \left(k\right)$ , $k=0,1,\cdots ,N-1,$ are added into the model.

Table 2 shows the simulation result obtained by using the algorithm proposed. It is mentioned that the minimum cost for the linear model-based optimal control problem is 5.9589 units. At the beginning of the iteration calculation procedure, the initial cost is 0.147463 unit, and a 90% of cost reduction is addressed to give the final cost of 0.014167 unit.

Number of iteration	Initial cost	Final cost	Elapsed time (s)
9	0.147463	0.014167	4.60934

Table 2.

Simulation result, Example 2.

The trajectories of the final control and the final state are, respectively, shown in Figures 7 and 8. It is noted that the state reaches to the steady state after 40 units of time by associating the control effort taken. This situation indicates that the temperature and the concentration are maintained at their steady state. Thus, the desired objective is confirmed. Figure 9 shows the costate behavior, which is reduced gradually to zero at the terminal time, and Figure 10 shows the stationary condition, which examines the existing of the optimal solution. The adjusted parameters, which are shown in Figures 11 and 12, respectively, measure the differences between the model used and the real plant.

Figure 7.

Final control u(k), Example 2.

Figure 8.

Final state x(k), Example 2.

Figure 9.

Final costate p(k), Example 2.

Figure 10.

Stationary H_u(k), Example 2.

Figure 11.

Adjusted parameter α(k), Example 2.

Figure 12.

Adjusted parameter γ(k), Example 2.

Hence, the correct optimal solution of Problem (P) is approximated successfully by solving the model in Problem (M), and the efficiency of the algorithm proposed is demonstrated.

5. Concluding remarks

The approach, which integrates system optimization and parameter estimation, was discussed in this chapter. The use of the conjugate gradient method in solving the model-based optimal control problem has been examined, and the applicability of the conjugate gradient approach in associating the principle of model-reality differences was identified. Definitely, many computational approaches could be used to solve the model-based optimal control; however, the algorithm proposed in this chapter gives a tractable solution procedure for handling the optimal control problems with different structures and parameters, especially for obtaining the optimal solution for the nonlinear optimal control problem. In conclusion, the efficiency of the algorithm is highly recommended. In future research, it is strongly suggested to investigate the application of optimization techniques in stochastic optimization and control.

Acknowledgments

The authors would like to acknowledge the Universiti Tun Hussein Onn Malaysia (UTHM) and the Ministry of Higher Education (MOHE) for the financial support for this study under the research grant FRGS VOT. 1561.

Conflict of interest

The authors declare no conflict of interest.