# Efficient Output Solution for Nonlinear Stochastic Optimal Control Problem with Model-Reality Differences

Authors: Sie Long Kek; Kok Lay Teo; Mohd Ismail Abdul Aziz
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A computational approach is proposed for solving the discrete time nonlinear stochastic optimal control problem. Our aim is to obtain the optimal output solution of the original optimal control problem through solving the simplified model-based optimal control problem iteratively. In our approach, the adjusted parameters are introduced into the model used such that the differences between the real system and the model used can be computed. Particularly, system optimization and parameter estimation are integrated interactively. On the other hand, the output is measured from the real plant and is fed back into the parameter estimation problem to establish a matching scheme. During the calculation procedure, the iterative solution is updated in order to approximate the true optimal solution of the original optimal control problem despite model-reality differences. For illustration, a wastewater treatment problem is studied and the results show the efficiency of the approach proposed.


---

## Body

## 1. Introduction

Many real world problems can be formulated as the stochastic dynamical systems [1–3]. In presence of the random noises, the exact state trajectory is impossible to be obtained. The output sequence, which is measured from the process plant, is also disturbed unavoidably. Since the fluctuation behavior of the output sequence would be the actual outcome, obtaining such outcome from a mathematical model is a challenging task. In stochastic system, estimation, identification, and adaptive control are the general techniques [4]. Particularly, the Kalman filtering theory and the extended Kalman filter give a great impact in studying stochastic systems, both for linear and nonlinear cases [5–8]. The data-driven method that could be applied in the fault diagnosis provides an efficient identification approach for stochastic systems. Using the process data to identify the parameters without knowing the actual process model is one of the advantages in modeling stochastic systems [9]. In addition, the stochastic switching systems, subject to random abrupt changes in their dynamics, attract the researchers to design, model, control, and optimize the stochastic systems [10, 11].The use of stochastic systems, therefore, plays the important role in the real world applications. The development of solution methods and the corresponding practical analysis are contributed to the stochastic research communities, ranging from engineering to business. From the literature, the applications of stochastic system have been well defined; see, for example, power management [12], portfolio selection [13], financial market debt crises [14], insurance with bankruptcy return [15], annuity contracts [16], natural gas networks [17], brain-machine interface operation [18], multi-degree-of-freedom systems [19], fleet composition problem [20], fault diagnosis [21], network control system [22], and stochastic switching system [23–25].In this paper, we propose a computational approach for optimal control of the nonlinear stochastic dynamical system in discrete time. Our aim is to obtain the optimal output solution of the original optimal control problem from a mathematical model. In doing so, a model-based optimal control problem is simplified from the original optimal control problem. Furthermore, the adjusted parameters are introduced into the model used. In this way, the differences between the real system and the model used can be computed. Thus, system optimization and parameter estimation are integrated interactively. On the other hand, the output, which is measured from the real plant, is fed back into the parameter estimation problem. This operation is implemented to establish a matching scheme, in turn, updating the optimal solution of the model used at each iteration step. Notice that the application of this operation, which is the advantage of the algorithm proposed, is in contrast to the works discussed in [26–28], where the real output is fed back into the system optimization problem. When convergence is achieved, the iterative solution approximates to the true optimal solution of the original optimal control problem, in spite of model-reality differences. Hence, the efficiency of the approach proposed is highly recommended.The rest of the paper is organized as follows. In Section2, a discrete time nonlinear stochastic optimal control problem is described and the corresponding model-based optimal control problem is simplified. In Section 3, an expanded optimal control model, which integrates system optimization and parameter estimation interactively, is introduced. Then, the iterative algorithm based on principle of model-reality differences is derived, and the computation procedure is summarized. In Section 4, a convergence analysis is provided. In Section 5, an example of the optimal control of a wastewater treatment problem is illustrated. Finally, some concluding remarks are made.


## 2. Problem Statement

Consider the following discrete time nonlinear stochastic optimal control problem:(1)
min
⁡
u
k
J
0
(
u
)
=
E
φ
x
N
,
N
+
∑
k
=
0
k
-
1
L
x
k
,
u
k
,
k
subject
to
x
k
+
1
=
f
x
k
,
u
k
,
k
+
G
ω
k
,
y
k
=
h
x
k
,
k
+
η
k
,
where u
(
k
)
∈
R
m, k
=
0,1
,
…
,
N
-
1, x
(
k
)
∈
R
n, k
=
0,1
,
…
,
N, and y
(
k
)
∈
R
p, k
=
0,1
,
…
,
N, are, respectively, the control sequence, the state sequence, and the output sequence. ω
(
k
)
∈
R
q, k
=
0,1
,
…
,
N
-
1, and η
(
k
)
∈
R
p, k
=
0,1
,
…
,
N, are the stationary Gaussian white noise sequences with zero mean and their covariances are given by Q
ω
∈
R
q
×
q and R
η
∈
R
p
×
p, which are positive definite matrices. G
∈
R
n
×
q is a process coefficient matrix, f
:
R
n
×
R
m
×
R
→
R
n represents the real plant, and h
:
R
n
×
R
→
R
p is the output measurement. J
0 is the scalar cost function and E
[
·
] is the expectation operator, whereas φ
:
R
n
×
R
→
R is the terminal cost and L
:
R
n
×
R
m
×
R
→
R is the cost under summation. It is assumed that all functions in (1) are continuously differentiable with respect to their respective arguments.The initial state is(2)
x
0
=
x
0
,
where x
0
∈
R
n is a random vector with mean and covariance given, respectively, by
(3)
E
[
x
0
]
=
x
-
0
,
E
x
0
-
x
-
0
x
0
-
x
-
0
T
=
M
0
.
Here, M
0
∈
R
n
×
n is a positive definite matrix. It is assumed that initial state, process noise, and measurement noise are statistically independent.This problem is referred to as Problem(
P
).Because of the complexity in the structure of the real plant and the presence of the random sequences, the exact solution of Problem(
P
) is impossible to be obtained. Moreover, applying the nonlinear filtering theory to estimate the state dynamics is computationally demanding. In view of these, we propose a simplified model-based optimal control problem, which is constructed by carrying out the linearization of Problem (
P
), in order to approximate the correct optimal solution of the original optimal control problem iteratively. This simplified model-based optimal control problem, which is referred to as Problem (
M
), is given by
(4)
min
⁡
u
k
J
1
u
=
1
2
x
-
N
T
S
N
x
-
N
+
γ
N
+
∑
k
=
0
N
-
1
1
2
x
-
k
T
Q
x
-
k
+
u
k
T
R
u
k
+
γ
k
subject
to
x
-
k
+
1
=
A
x
-
k
+
B
u
k
+
α
1
k
,
y
-
k
=
C
x
-
k
+
α
2
k
,
x
-
(
0
)
=
x
-
0
,
where x
-
(
k
)
∈
R
n
,
k
=
0,1
,
…
,
N, and y
-
(
k
)
∈
R
p
,
k
=
0,1
,
…
,
N, are the expected state sequence and the expected output sequence, respectively. α
1
(
k
)
∈
R
n
,
k
=
0,1
,
…
,
N
-
1, α
2
(
k
)
∈
R
p
,
k
=
0,1
,
…
,
N, and γ
(
k
)
∈
R
,
k
=
0,1
,
…
,
N, are adjustable parameters. A
∈
R
n
×
n is a state transition matrix, B
∈
R
n
×
m is a control coefficient matrix, and C
∈
R
p
×
n is an output coefficient matrix, while S
(
N
)
∈
R
n
×
n and Q
∈
R
n
×
n are positive semidefinite matrices and R
∈
R
m
×
m is a positive definite matrix.Notice that, without the adjustable parameters, Problem(
M
) is a standard linear quadratic regulator (LQR) optimal control problem. Solving this problem will not give us the optimal solution of the original optimal control problem. However, by adding the adjustable parameters into the model used, the differences between the real system and the model used can be computed such that system optimization and parameter estimation are integrated interactively.


## 3. System Optimization with Parameter Estimation

Now, let us introduce an expanded optimal control problem, which is referred to as Problem(
E
), given as follows:
(5)
min
⁡
u
k
J
2
(
u
)
=
1
2
x
-
N
T
S
N
x
-
N
+
γ
N
+
∑
k
=
0
N
-
1
1
2
x
-
k
T
Q
x
-
k
+
u
k
T
R
u
k
+
γ
k
+
1
2
r
1
u
k
-
v
k
2
+
1
2
r
2
x
-
k
-
z
k
2
subject
to
x
-
k
+
1
=
A
x
-
k
+
B
u
k
+
α
1
k
,
y
-
k
=
C
x
-
k
+
α
2
k
,
x
-
0
=
x
-
0
,
1
2
z
N
T
S
N
z
N
+
γ
N
=
φ
z
N
,
N
,
1
2
z
k
T
Q
z
k
+
v
k
T
R
v
k
+
γ
k
=
L
z
k
,
v
k
,
k
,
A
z
k
+
B
v
k
+
α
1
k
=
f
z
k
,
v
k
,
k
,
C
z
k
+
α
2
k
=
y
k
,
v
k
=
u
k
,
z
k
=
x
-
k
,
where v
(
k
)
∈
R
m
,
k
=
0,1
,
…
,
N
-
1, and z
(
k
)
∈
R
n
,
k
=
0,1
,
…
,
N, are introduced to separate the control and the expected state from the respective signals in the parameter estimation problem and · denotes the usual Euclidean norm. The terms 1
/
2
r
1
u
k
-
v
k
2 and 1
/
2
r
2
x
-
k
-
z
k
2 are introduced to improve convexity and enhance convergence of the resulting iterative algorithm. It is important to note that the algorithm is to be designed such that the constraints v
(
k
)
=
u
(
k
) and z
(
k
)
=
x
-
(
k
) will be satisfied at the end of the iterations. In this situation, the state estimate z
(
k
) and the control v
(
k
) will be used for the computation in the parameter estimation and the matching schemes. On the other hand, the corresponding expected state x
-
(
k
) and control u
(
k
) will be reserved for optimizing the model-based optimal control problem.It is important to note that the output measured from the real plant is fed back into the parameter estimation problem and the matching scheme, which aims at updating the model output from the model-based optimal control problem repeatedly. On this basis, the output residual could be reduced such that the model output approximates closely to the real output, in spite of model-reality differences. This improvement enhances the accuracy of the output solution as discussed in [26–28].


### 3.1. Optimality Conditions

Define the Hamiltonian function for Problem(
E
) as follows [28–30]:
(6)
H
k
=
1
2
x
-
k
T
Q
x
-
k
+
u
k
T
R
u
k
+
γ
k
+
1
2
r
1
u
k
-
v
k
2
+
1
2
r
2
x
-
k
-
z
k
2
-
λ
k
T
u
k
-
β
k
T
x
-
k
+
q
k
T
C
x
-
k
+
α
2
k
-
y
-
k
+
p
k
+
1
T
A
x
-
k
+
B
u
k
+
α
1
k
.
Then, the augmented cost function becomes
(7)
J
2
′
k
=
1
2
x
-
N
T
S
N
x
-
N
+
γ
(
N
)
+
p
0
T
x
-
0
-
p
N
T
x
-
N
+
ξ
N
×
φ
z
N
,
N
-
1
2
z
N
T
S
N
z
N
-
γ
N
+
Γ
T
x
-
N
-
z
N
+

∑
k
=
0
N
-
1
H
k
-
p
k
T
x
-
k

+
λ
k
T
v
k

+
β
k
T
z
k
+
ξ
k
1
2
L
z
k
,
v
k
,
k
-
1
2
z
k
T
Q
z
k
+
v
k
T
R
v
k
-
γ
k

1
2
+
μ
k
T
f
z
k
,
v
k
,
k
-
A
z
k
-
B
v
k
-
α
1
k
+
π
k
T
y
k
-
C
z
k
-
α
2
k
,
where p
(
k
)
,
q
(
k
)
,
ξ
(
k
)
,
μ
(
k
)
,
π
(
k
)
,
Γ
,
λ
(
k
), and β
(
k
) are the appropriate multipliers to be determined later.Applying the calculus of variation [26, 27, 29–31] to (7), the following necessary optimality conditions are obtained.(a)
Stationary condition:(8a)
R
u
k
+
B
T
p
k
+
1
-
λ
k
-
r
1
v
k
-
u
k
=
0
.
(b)
Costate equation:(8b)
p
k
=
Q
x
-
k
+
A
T
p
k
+
1
-
β
k
-
r
2
z
k
-
x
-
k
.
(c)
State equation:(8c)
x
-
(
k
+
1
)
=
A
x
-
(
k
)
+
B
u
(
k
)
+
α
1
(
k
)with the boundary conditions x
-
(
0
)
=
x
-
0 and p
(
N
)
=
Γ.
(d)
Adjustable parameter equations:(9a)
φ
z
N
,
N
=
1
2
z
N
T
S
N
z
N
+
γ
N
,
(9b)
L
z
k
,
v
k
,
k
=
1
2
z
k
T
Q
z
k
+
v
k
T
R
v
k
+
γ
k
,
(9c)
f
z
k
,
v
k
,
k
=
A
z
k
+
B
v
k
+
α
1
k
,
(9d)
y
k
=
C
z
k
+
α
2
k
.
(e)
Multiplier equations:(10a)
Γ
=
∇
z
(
k
)
φ
-
S
N
z
N
,
(10b)
λ
k
=
-
∇
v
k
L
-
R
v
k
-
∂
f
∂
v
k
-
B
T
p
^
k
+
1
,
(10c)
β
k
=
-
∇
z
k
L
-
Q
z
k
-
∂
f
∂
z
k
-
A
T
p
^
k
+
1with ξ
(
k
)
=
1
,
μ
(
k
)
=
p
^
(
k
+
1
), and π
(
k
)
=
q
(
k
)
=
0.
(f)
Separable variables:(11)
v
k
=
u
k
,
z
k
=
x
-
k
,
p
^
k
=
p
k
.In view of these necessary optimality conditions, conditions (8a), (8b), and (8c) are the necessary conditions for the modified model-based optimal control problem, conditions (9a), (9b), (9c), and (9d) define the parameter estimation problem, and conditions (10a), (10b), and (10c) are used to compute the multipliers.


### 3.2. Feedback Control Law

Taking the necessary optimality conditions (8a), (8b), and (8c), the modified model-based optimal control problem, which is referred to as Problem (
M
M
), is defined as follows:
(12)
min
⁡
u
k
J
3
(
u
)
=
1
2
x
-
N
T
S
(
N
)
x
-
(
N
)
+
γ
(
N
)
+
Γ
T
x
-
N
+
∑
k
=
0
N
-
1
1
2
x
-
k
T
Q
x
-
k
+
u
k
T
R
u
k
+
γ
(
k
)
+
1
2
r
1
u
k
-
v
k
2
+
1
2
r
2
x
-
k
-
z
k
2
-
λ
k
T
u
k
-
β
k
T
x
-
k
subject
to
x
-
k
+
1
=
A
x
-
k
+
B
u
k
+
α
1
k
,
y
-
k
=
C
x
-
k
+
α
2
k
,
x
-
(
0
)
=
x
-
0
.To solve Problem(
M
M
), we will construct a feedback control law, which includes the model-reality differences for system optimization. Hence, with the determined value of the adjustable parameters, the corresponding result is stated in the following theorem.Theorem 1 (expanded optimal control law).
Suppose that the expanded optimal control law for Problem(
E
) exists. Then, this optimal control law is the feedback control law for Problem (
M
M
) given by
(13)
u
k
=
-
K
k
x
-
k
+
u
f
f
k
,
where(14a)
u
f
f
k
=
-
R
a
+
B
T
S
k
+
1
B
-
1
×
B
T
s
k
+
1
+
B
T
S
k
+
1
α
1
k
-
λ
a
k
,
(14b)
K
k
=
R
a
+
B
T
S
k
+
1
B
-
1
B
T
S
k
+
1
A
,
(14c)
S
k
=
Q
a
+
A
T
S
k
+
1
A
-
B
K
k
,
(14d)
s
k
=
A
-
B
K
k
T
s
k
+
1
+
S
k
+
1
α
1
k
+
K
k
T
λ
a
k
-
β
a
kwith the boundary conditions S
(
N
) given and s
(
N
)
=
0, and
(15)
R
a
=
R
+
r
1
I
m
;
Q
a
=
Q
+
r
2
I
n
;
λ
a
(
k
)
=
λ
(
k
)
+
r
1
v
(
k
)
;
β
a
(
k
)
=
β
(
k
)
+
r
2
z
(
k
)
.Proof.
From (8a), the stationary condition can be rearranged by
(16)
R
a
u
k
=
-
B
T
p
k
+
1
+
λ
a
k
.
Applying the sweep method [28, 30, 31], that is,
(17)
p
k
=
S
k
x
-
k
+
s
k
,
substitute (17) for k
=
k
+
1 into (16) to yield
(18)
R
a
u
k
=
-
B
T
S
k
+
1
x
-
k
+
1
-
B
T
s
k
+
1
+
λ
a
k
.
Then, consider the state equation (8c) in (18). After some algebraic manipulations, the feedback control law (13) is obtained, where (14a) and (14b) are satisfied.
From (8b), the costate equation is rewritten as follows after substituting (17) for k
=
k
+
1 into (8b):
(19)
p
k
=
Q
a
x
-
k
+
A
T
S
k
+
1
x
-
k
+
1
+
A
T
s
k
+
1
-
β
a
k
.
Considering the state equation (8c) in (19), we have
(20)
p
k
=
Q
a
x
-
(
k
)
+
A
T
S
k
+
1
A
x
-
k
+
B
u
k
+
α
1
k
+
A
T
s
k
+
1
-
β
a
k
.
Use the feedback control law (13) in (20), and, doing some algebraic manipulations by considering (14a) and (14b), it is found that (14c) and (14d) are satisfied after comparing the manipulation result to (17). This completes the proof.Taking (13) in (8c), the state equation becomes
(21)
x
-
(
k
+
1
)
=
A
-
B
K
k
x
-
(
k
)
+
B
u
f
f
k
+
α
1
k
and the output is measured from
(22)
y
-
k
=
C
x
-
k
+
α
2
k
.


### 3.3. Adjustable Parameters and Multipliers

Now, we apply the separable variables given in (11) for solving the parameter estimation problem as defined in (9a), (9b), (9c), and (9d). Our aim is to establish the matching scheme, where the differences between the real system and the model used are taken into account. Consequently, the adjusted parameters, which are resulting from parameter estimation problem defined in (9a), (9b), (9c), and (9d), are calculated from(23a)
α
1
k
=
f
z
k
,
v
k
,
k
-
A
z
k
-
B
v
k
,
(23b)
α
2
k
=
y
k
-
C
z
k
,
(23c)
γ
N
=
φ
z
N
,
N
-
1
2
z
N
T
S
N
z
N
,
(23d)
γ
k
=
L
z
k
,
v
k
,
k
-
1
2
z
k
T
Q
z
k
+
v
k
T
R
v
k
.The multipliers, which are related to the Jacobian matrix of the functionsf and L with respect to v
(
k
) and z
(
k
), are computed from (24a)
Γ
=
∇
z
(
k
)
φ
-
S
N
z
N
,
(24b)
λ
k
=
-
∇
v
k
L
-
R
v
k
-
∂
f
∂
v
k
-
B
T
p
^
k
+
1
,
(24c)
β
k
=
-
∇
z
k
L
-
Q
z
k
-
∂
f
∂
z
k
-
A
T
p
^
k
+
1
.


### 3.4. Iterative Algorithm

From the discussion above, the resulting algorithm, which is an iterative algorithm, is summarized below.Iterative AlgorithmData. Consider A, B, C, G, Q, R, Q
ω, R
η, S
(
N
), M
0, x
-
0, N, r
1, r
2, k
v, k
z, k
p, f, L, h, and φ. Note that A and B may be chosen through the linearization of f, and C is obtained from the linearization of h.Step 0. Compute a nominal solution. Assume α
1
(
k
)
=
0, k
=
0,1
,
…
,
N
-
1, α
2
(
k
)
=
0, γ
(
k
)
=
0, k
=
0,1
,
…
,
N, and r
1
=
r
2
=
0. Solve Problem (
M
) defined by (4) to obtain u
(
k
)
0, k
=
0,1
,
…
,
N
-
1, and x
-
(
k
)
0
,
y
-
(
k
)
0
,
p
(
k
)
0, k
=
0,1
,
…
,
N. Then, with α
1
(
k
)
=
0, k
=
0,1
,
…
,
N
-
1, α
2
(
k
)
=
0, γ
(
k
)
=
0
,
k
=
0,1
,
…
,
N, and r
1, r
2 from data, compute K
(
k
) and S
(
k
) from (14b) and (14c), respectively. Set i
=
0, v
(
k
)
0
=
u
(
k
)
0, z
(
k
)
0
=
x
-
(
k
)
0, and p
^
(
k
)
0
=
p
(
k
)
0.Step 1. Compute the adjustable parameter α
1
(
k
)
i, k
=
0,1
,
…
,
N
-
1, α
2
(
k
)
i, and γ
(
k
)
i, k
=
0,1
,
…
,
N, from (23a), (23b), (23c), and (23d). This is called the parameter estimation step.Step 2. Compute the modifiers Γ
i, λ
(
k
)
i, and β
(
k
)
i, k
=
0,1
,
…
,
N
-
1, from (24a), (24b), and (24c). This requires the partial derivatives of f, h, and L with respect to v
(
k
)
i and z
(
k
)
i.Step 3. With the determined α
1
(
k
)
i, α
2
(
k
)
i, γ
(
k
)
i, Γ
i, λ
(
k
)
i, β
(
k
)
i, v
(
k
)
i, and z
(
k
)
i, solve Problem (
M
M
) defined by (12) using the result as given in Theorem 1. This is called the system optimization step.(3.1)
Solve (14d) backward to obtain s
(
k
)
i
,
k
=
0,1
,
…
,
N, and solve (14a), either backward or forward, to obtain u
f
f
(
k
)
i
,
k
=
0,1
,
…
,
N
-
1.
(3.2)
Use (13) to obtain the new control u
(
k
)
i
,
k
=
0,1
,
…
,
N
-
1.
(3.3)
Use (21) to obtain the new state x
-
(
k
)
i
,
k
=
0,1
,
…
,
N.
(3.4)
Use (17) to obtain the new costate p
(
k
)
i
,
k
=
0,1
,
…
,
N.
(3.5)
Use (22) to obtain the new output y
-
(
k
)
i
,
k
=
0,1
,
…
,
N.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:(25a)
v
k
i
+
1
=
v
k
i
+
k
v
u
k
i
-
v
k
i
,
(25b)
z
k
i
+
1
=
z
k
i
+
k
z
x
-
k
i
-
z
k
i
,
(25c)
p
^
k
i
+
1
=
p
^
k
i
+
k
p
p
k
i
-
p
^
k
i
,where k
v
,
k
z
,
k
p
∈
(
0,1
] are scalar gains. If v
(
k
)
i
+
1
=
v
(
k
)
i, k
=
0,1
,
…
,
N
-
1, and z
(
k
)
i
+
1
=
z
(
k
)
i, k
=
0,1
,
…
,
N, within a given tolerance, stop; else set i
=
i
+
1 and repeat from Steps 1–4.Remarks
(a)
The offline computation is done, as stated in Step 0, to calculateK
(
k
), k
=
0,1
,
…
,
N
-
1, and S
(
k
)
,
k
=
0,1
,
…
,
N, for the control law design. Then, these parameters are used for solving Problem (
M
) in Step 0 and for solving Problem (
M
M
) in Step 3, respectively.
(b)
The variablesα
1
(
k
)
i, α
2
(
k
)
i, γ
(
k
)
i, Γ
i, λ
(
k
)
i, β
(
k
)
i, and s
(
k
)
i are zero in Step 0. Their calculated values, α
1
(
k
)
i, α
2
(
k
)
i, and γ
(
k
)
iin Step 1, Γ
i, λ
(
k
)
i, and β
(
k
)
i in Step 2, and s
(
k
)
i in Step 3, change from iteration to iteration.
(c)
The driving inputu
f
f
(
k
) in (14a) corrects the differences between the real plant and the model used, and it also drives the controller given in (13).
(d)
Problem(
P
) does not need to be linear or to have a quadratic cost function.
(e)
The conditionsv
(
k
)
i
=
v
(
k
)
i
+
1 and z
(
k
)
i
=
z
(
k
)
i
+
1 are required to be satisfied for the converged optimal control sequence and the converged state estimate sequence. The following averaged 2 norms are computed and then they are compared with a given tolerance to verify the convergence of v
(
k
) and z
(
k
):(26a)
v
i
+
1
-
v
i
2
=
1
N
-
1
∑
k
=
0
N
-
1
v
k
i
+
1
-
v
k
i
1
/
2
,
(26b)
z
i
+
1
-
z
i
2
=
1
N
∑
k
=
0
N
z
k
i
+
1
-
z
k
i
1
/
2
.
(f)
The relaxation scalars(
k
v
,
k
z
,
k
p
) are the step sizes in the regulating convergence mechanism. They could be normally chosen as a certain value in (0, 1] but this choice may not result in the optimal number of iterations. It is important to note that the optimal choice of k
v
,
k
z
,
k
p
∈
(
0,1
] is problem dependent. It is required to run the algorithm (from Step 1 to Step 4) several times. They are initially chosen as k
v
=
k
z
=
k
p
=
1 for first run of the algorithm (Steps 1–4), and then the algorithm is run with different values chosen from 0.1 to 0.9. The value which provides the optimal number of iterations can then be determined. The parameters r
1 and r
2 are to enhance the convexity so as to improve the convergence of the algorithm.


## 3.1. Optimality Conditions

Define the Hamiltonian function for Problem(
E
) as follows [28–30]:
(6)
H
k
=
1
2
x
-
k
T
Q
x
-
k
+
u
k
T
R
u
k
+
γ
k
+
1
2
r
1
u
k
-
v
k
2
+
1
2
r
2
x
-
k
-
z
k
2
-
λ
k
T
u
k
-
β
k
T
x
-
k
+
q
k
T
C
x
-
k
+
α
2
k
-
y
-
k
+
p
k
+
1
T
A
x
-
k
+
B
u
k
+
α
1
k
.
Then, the augmented cost function becomes
(7)
J
2
′
k
=
1
2
x
-
N
T
S
N
x
-
N
+
γ
(
N
)
+
p
0
T
x
-
0
-
p
N
T
x
-
N
+
ξ
N
×
φ
z
N
,
N
-
1
2
z
N
T
S
N
z
N
-
γ
N
+
Γ
T
x
-
N
-
z
N
+

∑
k
=
0
N
-
1
H
k
-
p
k
T
x
-
k

+
λ
k
T
v
k

+
β
k
T
z
k
+
ξ
k
1
2
L
z
k
,
v
k
,
k
-
1
2
z
k
T
Q
z
k
+
v
k
T
R
v
k
-
γ
k

1
2
+
μ
k
T
f
z
k
,
v
k
,
k
-
A
z
k
-
B
v
k
-
α
1
k
+
π
k
T
y
k
-
C
z
k
-
α
2
k
,
where p
(
k
)
,
q
(
k
)
,
ξ
(
k
)
,
μ
(
k
)
,
π
(
k
)
,
Γ
,
λ
(
k
), and β
(
k
) are the appropriate multipliers to be determined later.Applying the calculus of variation [26, 27, 29–31] to (7), the following necessary optimality conditions are obtained.(a)
Stationary condition:(8a)
R
u
k
+
B
T
p
k
+
1
-
λ
k
-
r
1
v
k
-
u
k
=
0
.
(b)
Costate equation:(8b)
p
k
=
Q
x
-
k
+
A
T
p
k
+
1
-
β
k
-
r
2
z
k
-
x
-
k
.
(c)
State equation:(8c)
x
-
(
k
+
1
)
=
A
x
-
(
k
)
+
B
u
(
k
)
+
α
1
(
k
)with the boundary conditions x
-
(
0
)
=
x
-
0 and p
(
N
)
=
Γ.
(d)
Adjustable parameter equations:(9a)
φ
z
N
,
N
=
1
2
z
N
T
S
N
z
N
+
γ
N
,
(9b)
L
z
k
,
v
k
,
k
=
1
2
z
k
T
Q
z
k
+
v
k
T
R
v
k
+
γ
k
,
(9c)
f
z
k
,
v
k
,
k
=
A
z
k
+
B
v
k
+
α
1
k
,
(9d)
y
k
=
C
z
k
+
α
2
k
.
(e)
Multiplier equations:(10a)
Γ
=
∇
z
(
k
)
φ
-
S
N
z
N
,
(10b)
λ
k
=
-
∇
v
k
L
-
R
v
k
-
∂
f
∂
v
k
-
B
T
p
^
k
+
1
,
(10c)
β
k
=
-
∇
z
k
L
-
Q
z
k
-
∂
f
∂
z
k
-
A
T
p
^
k
+
1with ξ
(
k
)
=
1
,
μ
(
k
)
=
p
^
(
k
+
1
), and π
(
k
)
=
q
(
k
)
=
0.
(f)
Separable variables:(11)
v
k
=
u
k
,
z
k
=
x
-
k
,
p
^
k
=
p
k
.In view of these necessary optimality conditions, conditions (8a), (8b), and (8c) are the necessary conditions for the modified model-based optimal control problem, conditions (9a), (9b), (9c), and (9d) define the parameter estimation problem, and conditions (10a), (10b), and (10c) are used to compute the multipliers.


## 3.2. Feedback Control Law

Taking the necessary optimality conditions (8a), (8b), and (8c), the modified model-based optimal control problem, which is referred to as Problem (
M
M
), is defined as follows:
(12)
min
⁡
u
k
J
3
(
u
)
=
1
2
x
-
N
T
S
(
N
)
x
-
(
N
)
+
γ
(
N
)
+
Γ
T
x
-
N
+
∑
k
=
0
N
-
1
1
2
x
-
k
T
Q
x
-
k
+
u
k
T
R
u
k
+
γ
(
k
)
+
1
2
r
1
u
k
-
v
k
2
+
1
2
r
2
x
-
k
-
z
k
2
-
λ
k
T
u
k
-
β
k
T
x
-
k
subject
to
x
-
k
+
1
=
A
x
-
k
+
B
u
k
+
α
1
k
,
y
-
k
=
C
x
-
k
+
α
2
k
,
x
-
(
0
)
=
x
-
0
.To solve Problem(
M
M
), we will construct a feedback control law, which includes the model-reality differences for system optimization. Hence, with the determined value of the adjustable parameters, the corresponding result is stated in the following theorem.Theorem 1 (expanded optimal control law).
Suppose that the expanded optimal control law for Problem(
E
) exists. Then, this optimal control law is the feedback control law for Problem (
M
M
) given by
(13)
u
k
=
-
K
k
x
-
k
+
u
f
f
k
,
where(14a)
u
f
f
k
=
-
R
a
+
B
T
S
k
+
1
B
-
1
×
B
T
s
k
+
1
+
B
T
S
k
+
1
α
1
k
-
λ
a
k
,
(14b)
K
k
=
R
a
+
B
T
S
k
+
1
B
-
1
B
T
S
k
+
1
A
,
(14c)
S
k
=
Q
a
+
A
T
S
k
+
1
A
-
B
K
k
,
(14d)
s
k
=
A
-
B
K
k
T
s
k
+
1
+
S
k
+
1
α
1
k
+
K
k
T
λ
a
k
-
β
a
kwith the boundary conditions S
(
N
) given and s
(
N
)
=
0, and
(15)
R
a
=
R
+
r
1
I
m
;
Q
a
=
Q
+
r
2
I
n
;
λ
a
(
k
)
=
λ
(
k
)
+
r
1
v
(
k
)
;
β
a
(
k
)
=
β
(
k
)
+
r
2
z
(
k
)
.Proof.
From (8a), the stationary condition can be rearranged by
(16)
R
a
u
k
=
-
B
T
p
k
+
1
+
λ
a
k
.
Applying the sweep method [28, 30, 31], that is,
(17)
p
k
=
S
k
x
-
k
+
s
k
,
substitute (17) for k
=
k
+
1 into (16) to yield
(18)
R
a
u
k
=
-
B
T
S
k
+
1
x
-
k
+
1
-
B
T
s
k
+
1
+
λ
a
k
.
Then, consider the state equation (8c) in (18). After some algebraic manipulations, the feedback control law (13) is obtained, where (14a) and (14b) are satisfied.
From (8b), the costate equation is rewritten as follows after substituting (17) for k
=
k
+
1 into (8b):
(19)
p
k
=
Q
a
x
-
k
+
A
T
S
k
+
1
x
-
k
+
1
+
A
T
s
k
+
1
-
β
a
k
.
Considering the state equation (8c) in (19), we have
(20)
p
k
=
Q
a
x
-
(
k
)
+
A
T
S
k
+
1
A
x
-
k
+
B
u
k
+
α
1
k
+
A
T
s
k
+
1
-
β
a
k
.
Use the feedback control law (13) in (20), and, doing some algebraic manipulations by considering (14a) and (14b), it is found that (14c) and (14d) are satisfied after comparing the manipulation result to (17). This completes the proof.Taking (13) in (8c), the state equation becomes
(21)
x
-
(
k
+
1
)
=
A
-
B
K
k
x
-
(
k
)
+
B
u
f
f
k
+
α
1
k
and the output is measured from
(22)
y
-
k
=
C
x
-
k
+
α
2
k
.


## 3.3. Adjustable Parameters and Multipliers

Now, we apply the separable variables given in (11) for solving the parameter estimation problem as defined in (9a), (9b), (9c), and (9d). Our aim is to establish the matching scheme, where the differences between the real system and the model used are taken into account. Consequently, the adjusted parameters, which are resulting from parameter estimation problem defined in (9a), (9b), (9c), and (9d), are calculated from(23a)
α
1
k
=
f
z
k
,
v
k
,
k
-
A
z
k
-
B
v
k
,
(23b)
α
2
k
=
y
k
-
C
z
k
,
(23c)
γ
N
=
φ
z
N
,
N
-
1
2
z
N
T
S
N
z
N
,
(23d)
γ
k
=
L
z
k
,
v
k
,
k
-
1
2
z
k
T
Q
z
k
+
v
k
T
R
v
k
.The multipliers, which are related to the Jacobian matrix of the functionsf and L with respect to v
(
k
) and z
(
k
), are computed from (24a)
Γ
=
∇
z
(
k
)
φ
-
S
N
z
N
,
(24b)
λ
k
=
-
∇
v
k
L
-
R
v
k
-
∂
f
∂
v
k
-
B
T
p
^
k
+
1
,
(24c)
β
k
=
-
∇
z
k
L
-
Q
z
k
-
∂
f
∂
z
k
-
A
T
p
^
k
+
1
.


## 3.4. Iterative Algorithm

From the discussion above, the resulting algorithm, which is an iterative algorithm, is summarized below.Iterative AlgorithmData. Consider A, B, C, G, Q, R, Q
ω, R
η, S
(
N
), M
0, x
-
0, N, r
1, r
2, k
v, k
z, k
p, f, L, h, and φ. Note that A and B may be chosen through the linearization of f, and C is obtained from the linearization of h.Step 0. Compute a nominal solution. Assume α
1
(
k
)
=
0, k
=
0,1
,
…
,
N
-
1, α
2
(
k
)
=
0, γ
(
k
)
=
0, k
=
0,1
,
…
,
N, and r
1
=
r
2
=
0. Solve Problem (
M
) defined by (4) to obtain u
(
k
)
0, k
=
0,1
,
…
,
N
-
1, and x
-
(
k
)
0
,
y
-
(
k
)
0
,
p
(
k
)
0, k
=
0,1
,
…
,
N. Then, with α
1
(
k
)
=
0, k
=
0,1
,
…
,
N
-
1, α
2
(
k
)
=
0, γ
(
k
)
=
0
,
k
=
0,1
,
…
,
N, and r
1, r
2 from data, compute K
(
k
) and S
(
k
) from (14b) and (14c), respectively. Set i
=
0, v
(
k
)
0
=
u
(
k
)
0, z
(
k
)
0
=
x
-
(
k
)
0, and p
^
(
k
)
0
=
p
(
k
)
0.Step 1. Compute the adjustable parameter α
1
(
k
)
i, k
=
0,1
,
…
,
N
-
1, α
2
(
k
)
i, and γ
(
k
)
i, k
=
0,1
,
…
,
N, from (23a), (23b), (23c), and (23d). This is called the parameter estimation step.Step 2. Compute the modifiers Γ
i, λ
(
k
)
i, and β
(
k
)
i, k
=
0,1
,
…
,
N
-
1, from (24a), (24b), and (24c). This requires the partial derivatives of f, h, and L with respect to v
(
k
)
i and z
(
k
)
i.Step 3. With the determined α
1
(
k
)
i, α
2
(
k
)
i, γ
(
k
)
i, Γ
i, λ
(
k
)
i, β
(
k
)
i, v
(
k
)
i, and z
(
k
)
i, solve Problem (
M
M
) defined by (12) using the result as given in Theorem 1. This is called the system optimization step.(3.1)
Solve (14d) backward to obtain s
(
k
)
i
,
k
=
0,1
,
…
,
N, and solve (14a), either backward or forward, to obtain u
f
f
(
k
)
i
,
k
=
0,1
,
…
,
N
-
1.
(3.2)
Use (13) to obtain the new control u
(
k
)
i
,
k
=
0,1
,
…
,
N
-
1.
(3.3)
Use (21) to obtain the new state x
-
(
k
)
i
,
k
=
0,1
,
…
,
N.
(3.4)
Use (17) to obtain the new costate p
(
k
)
i
,
k
=
0,1
,
…
,
N.
(3.5)
Use (22) to obtain the new output y
-
(
k
)
i
,
k
=
0,1
,
…
,
N.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:(25a)
v
k
i
+
1
=
v
k
i
+
k
v
u
k
i
-
v
k
i
,
(25b)
z
k
i
+
1
=
z
k
i
+
k
z
x
-
k
i
-
z
k
i
,
(25c)
p
^
k
i
+
1
=
p
^
k
i
+
k
p
p
k
i
-
p
^
k
i
,where k
v
,
k
z
,
k
p
∈
(
0,1
] are scalar gains. If v
(
k
)
i
+
1
=
v
(
k
)
i, k
=
0,1
,
…
,
N
-
1, and z
(
k
)
i
+
1
=
z
(
k
)
i, k
=
0,1
,
…
,
N, within a given tolerance, stop; else set i
=
i
+
1 and repeat from Steps 1–4.Remarks
(a)
The offline computation is done, as stated in Step 0, to calculateK
(
k
), k
=
0,1
,
…
,
N
-
1, and S
(
k
)
,
k
=
0,1
,
…
,
N, for the control law design. Then, these parameters are used for solving Problem (
M
) in Step 0 and for solving Problem (
M
M
) in Step 3, respectively.
(b)
The variablesα
1
(
k
)
i, α
2
(
k
)
i, γ
(
k
)
i, Γ
i, λ
(
k
)
i, β
(
k
)
i, and s
(
k
)
i are zero in Step 0. Their calculated values, α
1
(
k
)
i, α
2
(
k
)
i, and γ
(
k
)
iin Step 1, Γ
i, λ
(
k
)
i, and β
(
k
)
i in Step 2, and s
(
k
)
i in Step 3, change from iteration to iteration.
(c)
The driving inputu
f
f
(
k
) in (14a) corrects the differences between the real plant and the model used, and it also drives the controller given in (13).
(d)
Problem(
P
) does not need to be linear or to have a quadratic cost function.
(e)
The conditionsv
(
k
)
i
=
v
(
k
)
i
+
1 and z
(
k
)
i
=
z
(
k
)
i
+
1 are required to be satisfied for the converged optimal control sequence and the converged state estimate sequence. The following averaged 2 norms are computed and then they are compared with a given tolerance to verify the convergence of v
(
k
) and z
(
k
):(26a)
v
i
+
1
-
v
i
2
=
1
N
-
1
∑
k
=
0
N
-
1
v
k
i
+
1
-
v
k
i
1
/
2
,
(26b)
z
i
+
1
-
z
i
2
=
1
N
∑
k
=
0
N
z
k
i
+
1
-
z
k
i
1
/
2
.
(f)
The relaxation scalars(
k
v
,
k
z
,
k
p
) are the step sizes in the regulating convergence mechanism. They could be normally chosen as a certain value in (0, 1] but this choice may not result in the optimal number of iterations. It is important to note that the optimal choice of k
v
,
k
z
,
k
p
∈
(
0,1
] is problem dependent. It is required to run the algorithm (from Step 1 to Step 4) several times. They are initially chosen as k
v
=
k
z
=
k
p
=
1 for first run of the algorithm (Steps 1–4), and then the algorithm is run with different values chosen from 0.1 to 0.9. The value which provides the optimal number of iterations can then be determined. The parameters r
1 and r
2 are to enhance the convexity so as to improve the convergence of the algorithm.


## 4. Convergence Analysis

In this section, we aim to provide a convergence analysis for the proposed computation procedure. The following assumptions are needed.Assumption 1.
(a) The derivatives off, L, and h exist; (b) the solution (
u
*
,
x
-
*
) is the exact optimal expected solution to Problem (
P
) with the optimal real output solution y
*.
The convergence property of the approximated solution to the true optimal solution is addressed in the following theorem.Theorem 2.
Let(
u
c
,
x
-
c
) be the converged solution to Problem (
M
). Then, there exist a converged output sequence y
-
c and a sequence of {
y
*
}, which is the original output sequence, such that
(27)
y
-
c
⟶
y
*
a
s
u
c
,
x
-
c
⟶
u
*
,
x
-
*
.
That is, the converged solution is the true optimal solution of Problem (
P
) and the converged output sequence is the true real output sequence.Proof.
Consider the real system of Problem(
P
) with the exact optimal expected solution (
u
*
,
x
-
*
) and the optimal real output solution y
* is given by (28a)
x
-
k
+
1
*
=
f
x
-
k
*
,
u
k
*
,
k
,
(28b)
y
-
k
*
=
h
x
-
k
*
,
k
,
(28c)
y
k
*
=
y
-
k
*
+
η
k
,where y
-
(
k
)
* is the exact optimal expected output solution and η
(
k
) is the output noise sequence. Meanwhile, in Problem (
M
), the model used consists of (29a)
x
-
k
+
1
c
=
A
x
-
k
c
+
B
u
k
c
+
α
1
k
i
,
(29b)
y
-
k
c
=
C
x
-
k
c
+
α
2
k
i
.Here, taking the adjusted parameters α
1
(
k
)
i and α
2
(
k
)
i from (23a) and (23b), the differences between the real system and the model used can be calculated from(30a)
α
1
k
i
=
f
z
k
i
,
v
k
i
,
k
-
A
z
k
i
-
B
v
k
i
,
(30b)
α
2
k
i
=
y
k
*
-
C
z
k
iat each iteration i. Note that (30a) and (30b) establish a matching scheme, in which, for any ρ
1
>
0, there exists a δ
1
>
0 such that, for k
=
0,1
,
…
,
N,
(31)
y
k
c
-
y
k
*
≤
ρ
1
whenever
u
c
,
x
-
c
-
v
i
,
z
i
<
δ
1
.
Hence, by substituting (30a) and (30b) into (29a) and (29b) and comparing the result yielded to (28a), (28b), and (28c), we conclude that
(32)
u
c
(
k
)
=
u
*
(
k
)
,
x
-
k
c
=
x
-
k
*
,
y
-
k
c
=
y
k
*
which are, respectively, the optimal expected solution and the optimal real output solution for the original optimal control problem. This completes the proof.


## 5. Illustrative Example

Consider the wastewater treatment problem [32–34]. The process equations, which are assumed to be unknown, are given by
(33)
x
1
k
+
1
x
2
k
+
1
=
1
+
μ
k
·
T
0
-
0.3636
·
T
1
x
1
k
x
2
k
+
-
T
·
x
1
k
0
-
T
·
x
2
k
T
u
1
k
u
1
(
k
)
·
u
2
k
+
ω
1
k
ω
2
k
y
(
k
)
=
x
1
(
k
)
+
η
(
k
)
with
(34)
μ
(
k
)
=
0.425
+
0.025
·
sin
(
2
π
k
/
96
)
·
x
2
(
k
)
0.4
+
x
2
(
k
)
,
for
k
=
0,1
,
…
,
95
,
where x
1
(
k
) is the methane gas flow rate, x
2
(
k
) is the substrate output concentration, u
1
(
k
) is the wastewater/dilution substance mix rate, u
2
(
k
) is the input flow rate, μ
(
k
) is the bacterial growth rate, and T is the sampling interval, which is 0.5 seconds. The initial state x
(
0
) has a mean given by x
-
0
=
(
0.5,1.6375
)
Tand its covariance is M
0
=
0
.
1
2
I
2, where I
2 is the two-dimensional identity matrix. The process noise ω
(
k
) and the measurement noise η
(
k
) have zero mean and their covariance is given by Q
ω
=
0
.
1
2
I
2 and R
η
=
0
.
1
2, respectively. Here, the aim is to determine an optimal control sequence u
(
k
)
=
(
u
1
(
k
)
,
u
2
(
k
)
)
T
∈
R
2 such that the cost function
(35)
J
0
u
=
1
2
E
x
1
288
2
+
x
2
288
2
+
∑
k
=
0
287
x
1
k
2
+
x
2
k
2
+
u
1
k
2
+
u
2
k
2

∑
k
=
0
287
is minimized subject to the dynamic system given by (34).This problem is regarded as Problem(
P
). The corresponding simplified model-based optimal control problem, which is referred to as Problem (
M
), is given by
(36)
min
⁡
u
k
J
1
(
u
)
=
∑
k
=
0
287
1
2
x
-
1
288
2
+
x
-
2
288
2
+
γ
288
+
∑
k
=
0
287
1
2
x
-
1
k
2
+
x
-
2
k
2
+
u
1
k
2
+
u
2
k
2
+
γ
k

∑
k
=
0
287
subject
to
x
-
(
k
+
1
)
=
1.0012
0
-
0.0018
0.9995
x
-
(
k
)
+
-
0.0025
0
-
0.0077
0.0005
u
k
+
α
1
k
y
-
(
k
)
=
x
-
1
(
k
)
+
α
2
(
k
)
with the initial state
(37)
x
-
1
(
0
)
=
0.5
x
-
2
(
0
)
=
1.6375
.By applying the algorithm proposed to solve Problem(
M
), the computation result is shown in Table 1. There is a 99.27 percent of reduction to the cost function, which gives the final cost 1.1810 units. The graphical results, which present the trajectories of output, state, and control, are shown, respectively, in Figures 1, 2, and 3. It is noticed that the model output sequence tracks closely to the real output sequence with the output residual 4.8
×
10
-
5. Both of the smooth trajectories of state and control show the optimal expected solution to the original optimal control problem.Table 1
Iteration result.
Iteration number
Elapsed time
Initial cost
Final cost
Output residual
15
2.708
160.8596
1.1810
0.000048Figure 1
Output trajectories: real outputy
(
k
) and model output y
-
(
k
).Figure 2
State trajectories: real statex
(
k
) and model state x
-
(
k
).Figure 3
Control trajectories:u
1
(
k
) and u
2
(
k
).Now, consider the target sequence that is a periodic square wave, with0.97,0.2
T for the first 48 time points and 1.0,0.1
T for the second time points, as discussed in [34]. Let this target sequence be the real state sequence in Problem (
P
). Then, the real output sequence is measured and is fed back into the parameter estimation problem. Here, the model used in Problem (
M
) remains the same as we mentioned above. Figure 4 shows the model output trajectory, which is generated by the algorithm proposed, and tracks the target sequence accordingly.Figure 4
Output trajectories: target sequencet
s
(
k
) and model output y
-
(
k
).From the results above, the output sequence, which is obtained by using the approach proposed, is efficient to optimal control of the discrete time nonlinear stochastic system. Hence, the applicability of the approach proposed is highlighted.


## 6. Concluding Remarks

In this paper, a computational approach was proposed, where the efficient output solution of the discrete time nonlinear stochastic optimal control problem is obtained. In our approach, the model-based optimal control problem is simplified from the original optimal control problem, where the adjusted parameters are introduced into the model used. On this basis, the differences between the real system and the model used are computed, and the integration of system optimization and parameter estimation could be made interactively. Establishment of the matching scheme by feeding back the output sequence that is measured from the real plant into the model used improves the accuracy of the model output sequence. For illustration, the wastewater treatment problem was studied, and the efficiency of the approach proposed is highly proven.



---

*Source: 659506-2015-03-02.xml*