- Email: [email protected]
A computational method for free terminal time optimal control problem governed by nonlinear time delayed systems
Accepted Manuscript
A computational method for free terminal time optimal control problem governed by nonlinear time delayed systems Qinqin Chai, Wu Wang PII: DOI: Reference:
S0307-904X(17)30541-3 10.1016/j.apm.2017.08.023 APM 11935
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
12 April 2016 18 August 2017 21 August 2017
Please cite this article as: Qinqin Chai, Wu Wang, A computational method for free terminal time optimal control problem governed by nonlinear time delayed systems, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.08.023
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • A free final time optimal control problem for nonlinear delay systems
CR IP T
is considered.
• A solving method based on fully informed particle swarm optimization is adopted.
AC
CE
PT
ED
M
AN US
• A practical optimal fishery control problem is solved.
1
ACCEPTED MANUSCRIPT
Qinqin Chaia,∗, Wu Wanga a
CR IP T
A computational method for free terminal time optimal control problem governed by nonlinear time delayed systems College of Electrical Engineering and Automation, Fuzhou University, Fuzhou, P. R. China
AN US
Abstract
This paper considers a free terminal time optimal control problem governed by nonlinear time delayed system, where both the terminal time and the control are required to be determined such that a cost function is minimized
M
subject to continuous inequality state constraints. To solve this free terminal time optimal control problem, the control parameterization technique is
ED
applied to approximate the control function as a piecewise constant control function, where both the heights and the switching times are regarded as
PT
decision variables. In this way, the free terminal time optimal control problem is approximated as a sequence of optimal parameter selection problems
CE
governed by nonlinear time delayed systems, each of which can be viewed as a nonlinear optimization problem. Then, a fully informed particle swarm optimization method is adopted to solve the approximate problem. Finally,
AC
two free terminal time optimal control problems, including an optimal fishery control problem, are solved by using the proposed method so as to demon∗
Corresponding author Email addresses: [email protected] (Qinqin Chai), [email protected] (Wu Wang) Preprint submitted to Applied Mathematical modelling
September 7, 2017
ACCEPTED MANUSCRIPT
strate its applicability. Keywords: time delayed system, free terminal time, optimal control
CR IP T
computation, fully informed particle swarm optimization, fishery control problem 1. Introduction
In this paper, we consider a class of free terminal time optimal control
AN US
problems governed by time delayed systems, where both the terminal times and control variables are to be chosen such that an objective function is minimized subject to some constraints. Many real world practical problems can be formulated in this form, see, for example, the resource (fishery, forestry) harvesting problem [1]. For fishery, the juvenile fish cannot be harvested,
M
and they need time to reach maturity for reproduction. This phenomenon is called time delay. Let the resource population be the state, and let the
ED
feeding, harvesting efforts be the controls. Then the growth of the resource population can be modelled as a dynamic system with delayed state [2]. Dif-
PT
ferent controls will result in different growth processes. Given the delayed dynamic system, the objective is to maximize a given economical benefit
CE
over a planning horizon, where the terminal time is the time when certain sustainable population is achieved. This unknown terminal time is called
AC
free terminal time. The corresponding optimal control problem, where the controls and the unknown terminal time are to be optimized simultaneously, is usually referred to as free terminal time optimal control problem. Other examples include (i) the biological process [3], where the feeding rates and the cells’ growth time are to be chosen such that, the yield is maximized and 3
ACCEPTED MANUSCRIPT
the operation cost is minimized; (ii) the HIV infection treatment process [4], where drugs are prescribed at crucial time points (to be determined) such
CR IP T
that the time from which the patient can live a normal life is the shortest; and (iii) attitude control of micro air vehicles [5], where the flight distance before running out of fuel or battery is maximized.
The free terminal time optimal control problem can, in principle, be solved
by using Pontryagin Maximum Principle [6]. However, as the terminal time is
AN US
determined by a stopping criterion, the resulting boundary value problem is
extremely difficult to solve.Consequently, different numerical computational approaches have been proposed in the literature for solving free terminal time optimal control problems. In [1], a two-stage algorithm is proposed for solving a free terminal time optimal control problem with time delay, where
M
Guinn’s transformation is used to transform the delayed dynamic system to a higher dimensional dynamic system without delay. Then, a gradient descent
ED
algorithm is used to solve the transformed problem. In [4], an adapted iterative forward backward sweep method integrated with a gradient algorithm is
PT
applied to solve a HIV infection problem, where the transversality condition for the terminal time is used. In [7, 8], necessary conditions of optimality
CE
and the transversality condition are obtained for an optimal control problem involving fractional differential equations by using Lagrange multiplier technique. On this basis, a shooting type of algorithm is utilized to search for the
AC
optimal final time and the optimal controls of the fractional optimal control problem. In [9], a hybrid time-scaling transformation and gradient-based algorithms are proposed to search for optimal initial states and terminal time for batch fermentation process. In [10, 11], computational methods based on
4
ACCEPTED MANUSCRIPT
control parameterization method are developed to solve free terminal time optimal control problems with or without time-delayed arguments. In [12],
CR IP T
an optimal switching control problem involving multiple time-delayed systems is formulated as a free terminal time control problem, where the cost
on changing control is part of the cost function. Then, a computational
method developed based on the control parameterization technique is presented. In [13], a two-stage optimization strategy is proposed, in which the
AN US
terminal time is adjusted in the outer stage, while a sequence of fixed time problems is solved in the inner stage. In [14], a simulation-based approxi-
mate dynamic programming method is proposed to solve a free-end optimal feedback control problem of fed-batch reactors. However, for most of these methods, they are developed for optimal control problems without delays
M
in their dynamical systems, and require the gradient of the cost function. These gradient-based algorithms are local search methods. The quality of
ED
the solution obtained depends critically on the initial guess of the decision variables for the optimization process. Furthermore, for the calculation of
PT
the gradient, it requires to solve the dynamic system forward in time and the corresponding co-state system backward in time. Furthermore, it is required
CE
to calculate the integration of the gradient of the Hamiltonian function over the planning horizon. Although all these tasks are to be carried out numerically, these tasks are known to be rather time consuming. In many real
AC
industrial applications, only input-output data are available [15]. For these situations, it is impossible to derive the gradient formulas. Thus, it is required to seek derivative-free optimization algorithms for finding the controls and the terminal time.
5
ACCEPTED MANUSCRIPT
In this paper, we propose to adopt a numerical computational approach based on particle swarm optimization (PSO) algorithm. PSO [16] is a derivative-
CR IP T
free and population-based stochastic search algorithm, which is easy to implement and has been found to be robust and fast in solving many practical
problems in sciences and engineering [17]. It is known that these PSO-based optimization algorithms can get trapped in local minima during the later stage of the searching process. Therefore, extensive research has been car-
AN US
ried out by many researchers over the years for overcoming this drawback. Improvements achieved are basically in the following three aspects: topology design [18], parameter adjustment [19], and hybrid algorithm with other optimization algorithms [20]. Many improved particle swarm optimization algorithms can enhance the performance of the original method at the expense
M
of increasing the complexity of these algorithms. For the topology design, it determines the social information sharing mechanism among populations.
ED
Subsequently, it changes the effectiveness of PSO. With this motivation, in this paper, a basic PSO with fully informed topology proposed in [21] is uti-
PT
lized. It is well recognized that this topology is simple and is more robust than other topologies.
CE
The rest of the paper is organised as follows. Problem formulation is given in Section 2. The problem approximation through the application of the control parameterization technique is detailed in Section 3. In Section 4,
AC
an effective optimization method based on a fully-informed particle swarm optimization method is presented and two numerical examples are solved. Finally, some concluding remarks are made in Section 5.
6
ACCEPTED MANUSCRIPT
2. Problem statement Consider the following time delayed system t ∈ [0, T ],
t ≤ 0,
x(t) = φ(t),
(1)
CR IP T
˙ x(t) = f (t, x(t), x(t − α), u(t)),
(2)
where T > 0 is an unknown terminal time; x(t) = [x1 (t), . . . , xn (t)]> ∈ Rn
is the state vector; x(t − α) ∈ Rn is the delayed state vector; α is a given
AN US
time delay; u(t) = [u1 (t), . . . , um (t)]> ∈ Rm is an unknown control function. f = [f1 , . . . , fn ]> and φ(t) = [φ1 , . . . , φn ]> are given functions.
A control function u = [u1 , . . . , um ]> is said to be an admissible control if the following boundedness constraints are satisfied i = 1, . . . , m,
t ∈ [0, T ],
(3)
M
bi ≤ ui (t) ≤ ci ,
where 0 ≤ bi < ci , i = 1, . . . , m, and bi , ci , i = 1, . . . , m, are, respectively,
ED
the lower and upper bounds of the ith control variable ui (t). Let U be the class of all admissible controls.
PT
Furthermore, we impose the following constraint: Tmin ≤ T ≤ Tmax
(4)
CE
where Tmin and Tmax are the lower and upper bounds of the terminal time, respectively.
AC
We assume that the following conditions are satisfied.
Assumption 1. The function f is continuously differentiable with respect to x and u for each t ∈ [0, T ], and piecewise continuous with respect to t for
each (x, u) ∈ Rn × Rm . Furthermore, the function φ is twice continuously differentiable. 7
ACCEPTED MANUSCRIPT
Assumption 2. There exists a real number L1 > 0 such that ˆ , u)| ≤ L1 (1 + |x| + |ˆ |f (t, x, x x| + |u|),
ˆ , u) ∈ R × Rn × Rn × Rm , (t, x, x
CR IP T
where | · | denotes the Euclidean norm.
On the basis of Assumptions 1 and 2, the dynamic system (1)-(2) admits a unique solution corresponding to each control u satisfying (3)[22]. Let
u ∈ U.
AN US
x(·|u) denote the solution of the dynamic system (1)-(2) corresponding to
We further require that the following state constraints are satisfied. Φl (t, x(t|u)) ≥ 0,
l = 1, . . . , r,
t ∈ [0, T ],
(5)
where the functions Φl , l = 1, . . . , r, are continuously differentiable with
M
respect to x for each t ∈ [0, T ] and piecewise continuous with respect to t for each x ∈ Rn . as follows:
ED
A free terminal time optimal control problem can now be formally stated
PT
Problem (P). Given system (1)-(2), find a terminal time T and an admis-
CE
sible control u(t) such that the performance function Z T J(T, u) = Φ0 (T, x(T |u)) + L0 (t, x(t|u), u)dt,
(6)
0
AC
is minimized subject to constraints (3)-(5). 3. Problem approximation In this section, we shall apply the control parameterization technique [23] to approximate Problem (P) by a sequences of nonlinear programming 8
ACCEPTED MANUSCRIPT
problems. The control parameterization technique approximates the control by a piecewise-constant function or piecewise smooth function with possible
CR IP T
discontinuities at the N − 1 pre-assigned partition time points, where the
heights of the piecewise-constant function are regarded as decision variables. Let uN i (t)
=
N X
σiN,k χ[τk−1 ,τk ) (t),
t ∈ [τk−1 , τk ),
k=1
AN US
Here, χ[τk−1 ,τk ) (t) is the indicator function defined by 1, if t ∈ [τk−1 , τk ), χ[τk−1 ,τk ) (t) = 0, elsewhere
i = 1, . . . , m.
(7)
(8)
and τk , k = 1, . . . , N , are the control switching times satisfying 0 = τ0 < τ1 < · · · < τN −1 < τN = T . For each k = 1, . . . , N and i = 1, . . . , m,
M
σiN,k is the control height of the ith piecewise constant control on the time
interval [τk−1 , τk ). Denote σ N = [(σ N,1 )> , . . . , (σ N,N )> ]> , where σ N,k =
ED
N,k > ] , k = 1, . . . , N . [σ1N,k , . . . , σm
Let ∆τk = τk − τk−1 , where the convention of ∆τ0 = 0 is adopted. It is
CE
PT
assumed that mink=1,...,N ∆τk > 0. Clearly, τk =
k X
∆τk .
i=1
In the classical control parameterization method, ∆τk is fixed. In this
AC
case, only the control heights are decision variables. Thus, to obtain an accurate optimal control policy, the number of N is required to be very large, resulting in a large number of decision variables. Thus, a non-uniform control vector parameterization is introduced in [24, 25], where the time durations ∆τk ,k = 1, ..., N , are also taken as decision variables. Let ∆τ N = 9
ACCEPTED MANUSCRIPT
[∆τ1 , . . . , ∆τN ]> . Then ∆τk and σiN,k , k = 1, . . . , N , i = 1, . . . , m, are decision variables. Let T N denote the set of all such ∆τ N and let ΞN be the set
CR IP T
containing all such σ N . Substituting (7) into system (1)-(2) gives the following system N X k=1
f (t, x(t), x(t − α), σ N,k )χIk (t),
x(t) = φ(t),
t ∈ Ik ,
t ≤ 0,
k = 1, . . . , N
(9)
(10)
P Pk−1 ∆τi , ki=1 ∆τi ). where Ik = [ i=1
AN US
˙ x(t) =
Note that, for each pair (σ N , ∆τ N ) ∈ ΞN × T N , the dynamic system
(9)-(10) admits a unique solution, which is denoted by x(·|σ N , ∆τ N ). For the constraints (4)-(5), they become
k=1
∆τk ≤ Tmax ,
M
Tmin ≤
N X
ED
bi ≤ σik ≤ ci , i = 1, . . . , m,
Φl (t, x(t|σ N , ∆τ N ) ≥ 0,
k = 1, . . . , N,
l = 1, . . . , r,
t ∈ [0, τN ].
(11) (12) (13)
PT
The new objective function is
CE
˜ N , ∆τ N ) = Φ0 (τN , x(τN |σ N , ∆τ N )) J(σ N Z τk X L0 (t, x(t|σ N , ∆τ N ), σ N,k )dt. + k=1
(14)
τk−1
AC
We may now state the approximate optimal control problem with contin-
uous state inequality constraints as follows: Problem (PN ). Given system (9)-(10), find a time duration vector ∆τ N
and a control parameter vector σ N such that the objective function (14) is minimized subject to constraints (11)-(13). 10
ACCEPTED MANUSCRIPT
Note that unless the dynamic system (9)-(10) is linear and the objective function (14) is linear or quadratic, they are non-convex. Thus, Problem (PN )
CR IP T
is non-convex. For example, we consider the optimal harvesting problem in [27]. Let us just discuss the situation when the partition number N = 1. The ˜ 1 , ∆τ1 )) three-dimensional drawing with respect to each pairs of (σ 1 , ∆τ1 , J(σ are shown in Figure 1, where the pairs of (σ 1 , ∆τ1 ) corresponding to the white
blocks are not feasible solutions. Thus, the feasible domain of Problem (PN )
AN US
is non-convex. In the following section, we will adopt an optimization strat-
egy developed based on a particle swarm optimization algorithm to solve Problem (PN ).
−5
M
−10 −15 −20
ED
The value of the cost function J
0
−25 1
Figure
1:
PT
u(t) 2
The
3
three-dimensional
4
6
drawing
t
2
with
respect
0
to
each
pairs
of
˜ 1 , ∆τ1 )) (σ , ∆τ1 , J(σ
AC
CE
1
4. Problem solving 4.1. A fully-informed particle swarm optimization Particle swarm optimization (PSO) is inspired by the social and cooperative behaviour of flocks of birds searching for food. A potential solution 11
ACCEPTED MANUSCRIPT
is called particle. A population of potential solutions is called swarm. The swarm is initialised with randomly generated particles, and the search for an
CR IP T
optimal result is carried out by updating the particles’ velocities and positions at each iteration. Recently, more and more efforts are devoted to the
design of methods for updating velocities and positions so as to avoid getting trapped in a local optimum or losing important information. In [21],
a fully informed particle swarm optimization (FIPSO) is proposed. In a
FIPSO are updated as follows. → − − → → − − v j+1 = χ(→ v ji + ϕ( P m − X ji )), i − → − → − = X ji + → v j+1 , X j+1 i i
AN US
D-dimensional search space, the velocities and positions of the particles in
i = 1, . . . , M,
j = 1, . . . , M axiter (15) (16)
i
M
where M is the population size; M axiter is the maximum iteration number; − → → − v j ∈ RD and X j ∈ RD are the velocity and position of the ith particle in the i
ED
population at the jth iteration, respectively; χ = 0.7298 is the constriction → − coefficient for the velocity; ϕ = 4.1 is the acceleration coefficient; P m refers
AC
CE
PT
to the comprehensive information learned from the ith particle’s neighbours. → − In this deterministic model, P m is calculated as → − Pm =
P
N→ − − W(k)→ ϕk Pk P , → − k∈N W(k) ϕ k → − → − ϕ k = U 0, ϕmax /|N | , k∈N
(17) (18)
→ − where N is the best among the neighbors of the particle; P k is the best position found by individual k, k ∈ N ; W(k) is the fitness difference between
the best position of the particle and the current individual, or it is a given
12
ACCEPTED MANUSCRIPT
constant. Furthermore, the velocity and position must satisfy i = 1, . . . , M,
j = 1, . . . , M axiter
(19)
i = 1, . . . , M,
j = 1, . . . , M axiter
(20)
CR IP T
− v min ≤ → v ji ≤ v max , − → X min ≤ X ji ≤ X max ,
where v min and v max are, respectively, the lower and upper bounds of the
velocity; X min and X max are, respectively, the lower and upper bounds of the position.
AN US
In the fully informed neighborhood, all neighbors are considered as sources of influence. Thus, the size and the topology of neighbors determine how diverse the influence will be. Experimental results in [21] show that FIPSO with the square topology, where the reference to the particle’s own index is removed from the neighborhood, converges faster than FIPSO with other
M
topologies. Furthermore, it has better chance of finding the best optimum.
ED
Thus, we will use FIPSO to solve the approximate optimal control problem. 4.2. A computational method based on FIPSO
PT
It is known that PSO can only solve problems with bound constraints, such as constraints (11) and (12). For optimization problems with complex
CE
constraints, such as continuous state inequality constraints, they cannot be
AC
solved directly by PSO. Thus, we introduce the following penalty function
13
ACCEPTED MANUSCRIPT
[26] to handle these complex constraints.
+
r X l=1
ρl
Z
τN
N Z X k=1
τk−1
L0 (t, x(t|σ N , ∆τ N ), σ N,k )dt
min{Φl (t, x(t|σ N , ∆τ N )), 0}2 dt
0
τk
CR IP T
J˜ρ (τN , σ N , ∆τ N ) = Φ0 (τN , x(τN |σ N , ∆τ N )) +
+ ρr+1 min{τN − Tmax , 0}2 + min{Tmin − τN , 0}2
(21)
AN US
where ρl , l = 1, . . . , r + 1, are given large penalty factors. The last three terms in (21) represent the constraint violations. Problem (PN ) is now approximated as:
Problem (PN,ρ ). Given system (9)-(10), find a time duration vector ∆τ N and a control vector σ N such that the objective function (21) is minimized
M
subject to bound constraint (12).
ED
Clearly, Problem (PN ) is a more accurate approximation of Problem (P ) when the partition number N increases. By using arguments similar to those given in [23], it can be shown that the approximate optimal cost will converge
PT
to the true optimal cost as N → ∞. However, in practice, we cannot let N = ∞. Furthermore, when N increases the decision parameters of the
CE
approximated problem will grow by m times. From our extensive numerical experiments, it suggests that the improvement in the cost value for N larger
AC
than N = 20 appears to be highly insignificant. Thus, we start with a small N , then Problem (PN,ρ ) is solved by FIPSO to obtain the optimal control.
We then increase the value of N and re-calculate the optimal control of the corresponding Problem (PN,ρ ). We repeat this process until the reduction in
the cost value is negligible. Details are listed as follows. 14
ACCEPTED MANUSCRIPT
Step 1: Initialize N . Set the swarm size M , the neighbour particles for each particle, and the maximum iteration or terminal conditions of FIPSO.
CR IP T
Step 2: Initialize the velocity and position of each particle in the swarm and its neighbours. Let the iteration p = 1. Denote the swarm by Swarmp .
Step 3: For each particle in the swarm, execute the following step: Solve the state differential equation (9)-(10) with the initial condition forward
AN US
in time from t = 0 to t = τN . Let the solution obtained be denoted as x(·|σ N , ∆τ N ). Substitute x(·|σ N , ∆τ N ) into the objective function (21) to obtain the cost value for each particle.
Step 4: Update each particle’s personal best position and the global best of the swarm. If the terminal conditions of FIPSO are satisfied, then go to
M
Step 6, else let p := p + 1 and go to Step 5.
ED
Step 5: Update the velocities and positions of the swarm according to (15)-(16) to obtain a new swarm Swarmp+1 . Go to Step 3.
PT
Step 6: Obtain the best control parameter vector σ N ∗ , time duration vector ∆τ N ∗ and the minimum cost value J˜ρ (τ N,∗ , σ N,∗ , ∆τ N,∗ ), and the best
CE
swarm P ∗M . If J˜ρ (τ N,∗ , σ N,∗ , ∆τ N,∗ ) − J˜ρ (τ N +1,∗ , σ N +1,∗ , ∆τ N +1,∗ ) ≤ ,
where = 10−6 is a required precision, then stop. Else let N := N + s,
AC
where s is some positive integer. Go to Step 2.
15
ACCEPTED MANUSCRIPT
5. Numerical results 5.1. Example 1
CR IP T
The first example involves optimal fishery harvesting, where the growth
of the population has limited resources (food, water, space). Suppose the growth of the population is limited by a carrying capacity K. Let x be the
current fish population and r be the intrinsic growth rate. The dynamics
of the population from generation to generation is governed by a logistic x(t) ]. K
Moreover, considering that the
AN US
model in the form of x(t) ˙ = rx(t)[1 −
juvenile fish cannot be harvested, and they needs time to reach maturity for reproduction, the population growth model should contain time delay from their birth to maturity. Hence, the following delayed logistic model with
M
harvesting for fish growth is considered [27]:
x(t) ˙ = 3x(t)[1 − 0.2x(t − 0.5)] − u(t), t ≤ 0,
(22) (23)
ED
x(t) = 2,
t ∈ [0, T ]
where u(t) is the harvesting effort, r = 3, K = 5.
PT
Clearly, the fish population increases, initially. However, when it gets closer to the carrying capacity, it becomes saturated forcing the rate of growth
CE
to decrease. To maximize sustainable yield and to achieve the commercial purpose of the fishery, it is necessary to harvest the population at appropriate
AC
time points with suitable harvest effort in a long run. Hence, the following control objective function are proposed : Z T Φ0 (T, u(t)) = 0.1T 2 + e−βt (CE x(t)−1 u(t)3 − pu(t))dt 0
16
(24)
ACCEPTED MANUSCRIPT
where β = 0.05 is the discount rate of the economic benefit influenced by the market; CE = 0.2 is the harvesting cost, p = 2 is the fish price. Over all, the
CR IP T
first term of equation (24) is used to shorten the total investment time. The second term is expenditure minus income. To minimize the second term is to maximum the economic benefit. In addition, to avoid overfishing and avoid exceeding the carrying capacity limitation, the following constraints must be satisfied: t ∈ [0, T ]
AN US
2 ≤x(t) ≤ 4,
(25)
0 ≤u(t) ≤ 3.5
(26)
0.1 ≤ T ≤ 100
(27)
The optimal fishery harvesting problem may now be described as follows.
M
Given the time delayed system (22)-(23), find a terminal time T and a harvesting strategy u such that the objective functional (24) is minimized,
ED
subject to constraints (25)-(27). By the approximation method described in Section 3, the free terminal
PT
time optimal control problem is approximated by the following optimization problem. Given the dynamic system
CE
N X
˙ x(t) =
AC
k=1
x(t) = 2,
{3x(t)[1 − 0.2x(t − 0.5)] − σ k }χIk ,
t ∈ Ik ,
k = 1, . . . , N (28)
t ≤ 0,
(29)
P Pk N where Ik = [ k−1 and a control i=1 ∆τi , i=1 ∆τi ), find a positive vector ∆τ 17
ACCEPTED MANUSCRIPT
parameter vector σ N such that the cost functional N
N Z X k=1
τk
τk−1
e−0.05t [0.2x(t)−1 (σ k )3 − 2σ k ]dt
CR IP T
J˜ρ (τN , σ , ∆τ ) = N
N N nX o2 X + 0.1 ∆τk + ρ1 (min{ ∆τk − 0.1, 0})2 k=1
k=1
+ ρ1 (min{100 − τN
0
k=1
∆τk , 0})2
[(min{x(t) − 2, 0})2 + (min{4 − x(t), 0})2 ]dt
AN US
+ ρ2
Z
N X
(30)
is minimized subject to the bound constraints (11) and (12), where ρ1 and ρ2 are penalty factors.
M
We solve this problem by writing a Matlab program on a personal computer with the following configuration: Intel Core i5-4590 3.3GHz CPU, 8G
ED
RAM, 64-bit Window 7 operating system. In order to test the performance of the proposed method, we choose N = 9 as given in [27]. The penalty factors are ρ1 = 1000 and ρ2 = 105 . The performances of the FIPSO are evaluated
PT
by comparing with that of the comprehensive learning particle swarm optimizer (CLPSO) proposed in [28], genetic algorithm (GA) given in [29]. For
CE
all three optimization algorithms, the swarm size or population size are 30, and the maximum iteration is 1000. According to reference [21], for FIPSO,
AC
the square topology is used for choosing four neighbours, and W(k) is the fitness difference between the best position of the particle and the current individual. For GA, the selection rate is 0.9, cross rate is 0.8, and mutation rate is 0.1. The other parameters of CLPSO are as defined in [28]. The optimal results are summarized in Table 1 and Figure 2-Figure 4. 18
ACCEPTED MANUSCRIPT
Table 1: Comparisons results of Example 1
FIPSO
CLPSO
GA
reference [14]
optimal cost value
-26.032
-26.025
-25.620
-25.97
optimal terminal time
12.1221
12.1487
12.1448
12.1565
CR IP T
optimization algorithm
AN US
Table 1 shows that FIPSO solves the optimal control problem with a shorter
terminal time (12.12 years), and a smaller cost value. Compared with [27], the harvest time can be shorted by 10 days. Figure 2 demonstrates that the convergence speed of FIPSO is faster than that of CLPSO and GA optimization algorithm. And the optimal state and control of FIPSO have a tendency
M
of becoming stable. Hence, FIPSO is superior to the other two optimization
are satisfied. 10
x 10
4
FIPSO
cost function value
CE AC
CLPSO
GA
−24.5
PT
8
ED
methods in solving this optimal control problem. In addition, all constraints
6
−25
4
−25.5 −26 490
2
495
500
0 100
200
300
400 500 600 iteration number
700
800
900
1000
Figure 2: The convergence trajectories of FIPSO, GA, and CLPSO.
19
ACCEPTED MANUSCRIPT
CR IP T
4
x(t)
3.5 3
CLPSO FIPSO GA
2.5
0
2
4
6
8
t
10
12
AN US
2
M
Figure 3: The optimal state trajectories obtained by FIPSO, GA, and CLPSO, respectively
ED
3.5
PT
u(t)
3
2.5
FIPSO CLPSO GA
1.5
0
2
4
6
AC
CE
2
t
8
10
12
Figure 4: The optimal piecewise-constant controls obtained by FIPSO, GA, and CLPSO, respectively
20
ACCEPTED MANUSCRIPT
Table 2: Comparisons results of Example 2
CLPSO
GA
reference [1]
CR IP T
optimization algorithm FIPSO optimal cost value
5.287
5.425
8.875
5.37
optimal terminal time
1.349
1.1818
1.005
1.446
AN US
5.2. Example 2
The second example is the optimization problem given in [1]
min J = subject to the dynamics
T
0
(x2 (t) + u2 (t))dt + 0.1T
M
u,T
Z
ED
x(t) ˙ = x(t) + 0.8x(t − 0.15) + 0.9u(t),
t ∈ [0, T ]
(31)
(32)
t ≤ 0,
(33)
−10 ≤ u(t) ≤ 10,
(34)
PT
x(t) = 1,
We solve this problem with control partition numbers N = 6, using
CE
FIPSO, CLPSO, GA algorithms with the same parameter settings for solving example 1. Comparison results are given in Table 2. We can also see that
AC
FIPSO is more effectiveness than the other optimization algorithms. From the results given in Table 2, we can see that the minimum cost has been reduced by 1.6% and the minimum time has been reduced by 7% compared with that of [1].
21
ACCEPTED MANUSCRIPT
6. Conclusion In this paper, a class of free terminal time optimal control problem in-
CR IP T
volving time delayed system is investigated. To solve this problem, we apply the control parameterization method to yield an approximate optimization
problem. Then a simple effective optimization method is developed, based on a fully informed particle swarm optimization method, to solve the approx-
imate problem. Numerical results show that the proposed method is more
AN US
effective than other existing solving methods with shorter terminal time and better cost function values. Acknowledgment
M
The authors thank the anonymous reviewers and Mr. Jianhua Liu for their valuable comments and suggestions to improve this paper. This work
ED
was supported by the Natural Science Foundation of Fujian Province, China (Grant No. 2016J05154), the National Natural Science Foundation of China
PT
(Grant No. 61673116, 61503131).
CE
References
[1] A. Boccia, P. Falugi, H. Maurer, R. B. Vinter. Free time optimal control problems with time delays. IEEE Conference on Decision and Control,
AC
2013:520-525.
[2] M. Liu, C. Bai. Optimal harvesting of a stochastic logistic model with time delay. Journal of Nonlinear Science, 2015, 25(2):277-289.
22
ACCEPTED MANUSCRIPT
[3] C. Liu, Z. Gong. Modelling and optimal control of a time-delayed switched system in fed-batch process. Journal of the Franklin Institute,
CR IP T
2014, 351(2): 840-856. [4] I. Elmouki. Free terminal time optimal control problem for the treatment of HIV infection. Jawra Journal of the American Water Resources Association, 2016, 7(1):148-161.
AN US
[5] J. J. Ewoud, Q. Chu, G. C. D. Croon. Adaptive incremental nonlinear
dynamic inversion for attitude control of micro air vehicles.Journal of Guidance Control and Dynamics, 2016, 39(3): 450-461. [6] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, E. F.
science, 1962.
M
Mishchenko, The mathematical theory of time optimal processes, Inter-
ED
[7] R. K. Biswas, S. Sen. Free final time fractional optimal control problems. Journal of the Franklin Institute, 2014, 351(2):941-951.
PT
[8] S. Pooseh, R. Almeida, D. F. M. Torres. Fractional Order Optimal Control Problems with Free Terminal Time. Journal of Industrial and Man-
CE
agement Optimization, 2014, 10(2): 363-381. [9] J. Yuan, C. Liu, X. Zhang, J. Xie, E. Feng, H. Yin, Z. Xiu. Optimal
AC
control of a batch fermentation process with nonlinear time-delay and free terminal time and cost sensitivity constraint. Journal of Process Control, 2016, 44: 41-52.
[10] R. Li, Y. Shi, Finite-time optimal consensus control for second-order 23
ACCEPTED MANUSCRIPT
multi-agent systems. Journal of Industrial & Management Optimization, 2014, 10(3):929-943.
CR IP T
[11] C. Jiang, Q. Lin, C. Yu, K. L. Teo, G. R. Duan. An exact penalty method
for free terminal time optimal control problem with continuous inequality constraints. Journal of Optimization Theory and Applications, 2012, 154(1):30-53.
AN US
[12] Z. Gong, C. Liu, Y. Wang. Optimal control of switched systems with multiple time-delays and a cost on changing control. Journal of Industrial and Management Optimization, 2017, 13, doi:10.3934/jimo.2017042.
[13] J. M. Jeong, S. J. Son, Time optimal control of semilinear control sys-
M
tems involving time delays. Journal of Optimization Theory and Appli-
ED
cations, 2015, 165(3):793-811.
[14] C. V. Peroni, N. S. Kaisare, J. H. Lee. Optimal control of a fed-batch bioreactor using simulation-based approximate dynamic programming.
CE
790.
PT
IEEE Transactions on Control Systems Technology, 2005, 13(5): 786-
[15] T. J. J. V. D. Boom, J. B. Klaassens, R. Meiland. Real-time time-
AC
optimal control for a nonlinear container crane using a neural network. 2007:79-84.
[16] R. Eberhat, J. Kennedy, A new optimizer using particle swarm theory. Proceedings of 6th International Symposium On Mico Machine and Human Science, 1995: 39-43. 24
ACCEPTED MANUSCRIPT
[17] C. Ou, W. Lin. Comparison between PSO and GA for parameters optimization of PID controller. IEEE International Conference on Mecha-
CR IP T
tronics and Automation. 2006: 2471-2475. [18] Q. Liu, W. Wei, H. Yuan, Z. H. Zhan, Y. Li. Topology selection for particle swarm optimization. Information Sciences,2016, 363: 154-173.
[19] A. B. Hashemi, M. R. Meybodi. A note on the learning automata based
puting, 2011, 11(1):689-705.
AN US
algorithms for adaptive parameter selection in PSO. Applied Soft Com-
[20] H. Garg. A hybrid PSO-GA algorithm for constrained optimization problems. Applied Mathematics & Computation, 2016, 274(11):292-305.
M
[21] R. Mendes, J. Kennedy, O. Neves, The fully informed particle swarm: simpler, maybe better. IEEE Transactions on evolutionary computation,
ED
2004, 8(3):204-210
[22] N. U. Ahmed, Dynamic systems and control with applications, Singa-
PT
pore: World Scientific, 2006. [23] K. L. Teo, C. J. Goh, K. H. Wong, A unified computational approach
CE
to optimal control problems, Essex: Longman Scientific and Technical,
AC
1991.
[24] Y. Lei, S. Li, Q. Zhang, X. Zhang, Free final time optimal control based on non-uniform control vector parameterization. Journal of Systems Science and Mathematical Sciences, 2012, 32(3): 277-287.
25
ACCEPTED MANUSCRIPT
[25] Q. Lin , R. Loxton, K. L. Teo, The control parameterization method for nonlinear optimal control: A survey. Journal of Industrial & Manage-
CR IP T
ment Optimization, 2017, 10(1):275-309. [26] C. J. Yu, K. L. Teo, L. S. Zhang, Y. Q. Bai, On a refinement of the con-
vergence analysis for the new exact penalty function method for continuous inequality constrained optimization problem. Journal of Industrial
AN US
Management and Opti- mization, 2012, 8(2):485-491.
[27] C. Liu, R. Loxton, K. L. Teo. A computational method for solving timedelay optimal control problems with free terminal time. Systems & Control Letters, 2014, 72: 53-60.
[28] J. J. Liang, A. K. Qin, P. N. Suganthan, S. Baskar, Comprehensive learn-
M
ing particle swarm optimizer for global optimization of multimodal func-
281-295.
ED
tions. IEEE Transactions on Evolutionary Computation, 2006, 10(3):
PT
[29] Goldberg D E. Genetic algorithm in search, optimization, and machine
AC
CE
learning. 1989, xiii(7): 2104C2116.
26
A computational method for free terminal time optimal control problem governed by nonlinear time delayed systems Qinqin Chai, Wu Wang PII: DOI: Reference:
S0307-904X(17)30541-3 10.1016/j.apm.2017.08.023 APM 11935
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
12 April 2016 18 August 2017 21 August 2017
Please cite this article as: Qinqin Chai, Wu Wang, A computational method for free terminal time optimal control problem governed by nonlinear time delayed systems, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.08.023
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • A free final time optimal control problem for nonlinear delay systems
CR IP T
is considered.
• A solving method based on fully informed particle swarm optimization is adopted.
AC
CE
PT
ED
M
AN US
• A practical optimal fishery control problem is solved.
1
ACCEPTED MANUSCRIPT
Qinqin Chaia,∗, Wu Wanga a
CR IP T
A computational method for free terminal time optimal control problem governed by nonlinear time delayed systems College of Electrical Engineering and Automation, Fuzhou University, Fuzhou, P. R. China
AN US
Abstract
This paper considers a free terminal time optimal control problem governed by nonlinear time delayed system, where both the terminal time and the control are required to be determined such that a cost function is minimized
M
subject to continuous inequality state constraints. To solve this free terminal time optimal control problem, the control parameterization technique is
ED
applied to approximate the control function as a piecewise constant control function, where both the heights and the switching times are regarded as
PT
decision variables. In this way, the free terminal time optimal control problem is approximated as a sequence of optimal parameter selection problems
CE
governed by nonlinear time delayed systems, each of which can be viewed as a nonlinear optimization problem. Then, a fully informed particle swarm optimization method is adopted to solve the approximate problem. Finally,
AC
two free terminal time optimal control problems, including an optimal fishery control problem, are solved by using the proposed method so as to demon∗
Corresponding author Email addresses: [email protected] (Qinqin Chai), [email protected] (Wu Wang) Preprint submitted to Applied Mathematical modelling
September 7, 2017
ACCEPTED MANUSCRIPT
strate its applicability. Keywords: time delayed system, free terminal time, optimal control
CR IP T
computation, fully informed particle swarm optimization, fishery control problem 1. Introduction
In this paper, we consider a class of free terminal time optimal control
AN US
problems governed by time delayed systems, where both the terminal times and control variables are to be chosen such that an objective function is minimized subject to some constraints. Many real world practical problems can be formulated in this form, see, for example, the resource (fishery, forestry) harvesting problem [1]. For fishery, the juvenile fish cannot be harvested,
M
and they need time to reach maturity for reproduction. This phenomenon is called time delay. Let the resource population be the state, and let the
ED
feeding, harvesting efforts be the controls. Then the growth of the resource population can be modelled as a dynamic system with delayed state [2]. Dif-
PT
ferent controls will result in different growth processes. Given the delayed dynamic system, the objective is to maximize a given economical benefit
CE
over a planning horizon, where the terminal time is the time when certain sustainable population is achieved. This unknown terminal time is called
AC
free terminal time. The corresponding optimal control problem, where the controls and the unknown terminal time are to be optimized simultaneously, is usually referred to as free terminal time optimal control problem. Other examples include (i) the biological process [3], where the feeding rates and the cells’ growth time are to be chosen such that, the yield is maximized and 3
ACCEPTED MANUSCRIPT
the operation cost is minimized; (ii) the HIV infection treatment process [4], where drugs are prescribed at crucial time points (to be determined) such
CR IP T
that the time from which the patient can live a normal life is the shortest; and (iii) attitude control of micro air vehicles [5], where the flight distance before running out of fuel or battery is maximized.
The free terminal time optimal control problem can, in principle, be solved
by using Pontryagin Maximum Principle [6]. However, as the terminal time is
AN US
determined by a stopping criterion, the resulting boundary value problem is
extremely difficult to solve.Consequently, different numerical computational approaches have been proposed in the literature for solving free terminal time optimal control problems. In [1], a two-stage algorithm is proposed for solving a free terminal time optimal control problem with time delay, where
M
Guinn’s transformation is used to transform the delayed dynamic system to a higher dimensional dynamic system without delay. Then, a gradient descent
ED
algorithm is used to solve the transformed problem. In [4], an adapted iterative forward backward sweep method integrated with a gradient algorithm is
PT
applied to solve a HIV infection problem, where the transversality condition for the terminal time is used. In [7, 8], necessary conditions of optimality
CE
and the transversality condition are obtained for an optimal control problem involving fractional differential equations by using Lagrange multiplier technique. On this basis, a shooting type of algorithm is utilized to search for the
AC
optimal final time and the optimal controls of the fractional optimal control problem. In [9], a hybrid time-scaling transformation and gradient-based algorithms are proposed to search for optimal initial states and terminal time for batch fermentation process. In [10, 11], computational methods based on
4
ACCEPTED MANUSCRIPT
control parameterization method are developed to solve free terminal time optimal control problems with or without time-delayed arguments. In [12],
CR IP T
an optimal switching control problem involving multiple time-delayed systems is formulated as a free terminal time control problem, where the cost
on changing control is part of the cost function. Then, a computational
method developed based on the control parameterization technique is presented. In [13], a two-stage optimization strategy is proposed, in which the
AN US
terminal time is adjusted in the outer stage, while a sequence of fixed time problems is solved in the inner stage. In [14], a simulation-based approxi-
mate dynamic programming method is proposed to solve a free-end optimal feedback control problem of fed-batch reactors. However, for most of these methods, they are developed for optimal control problems without delays
M
in their dynamical systems, and require the gradient of the cost function. These gradient-based algorithms are local search methods. The quality of
ED
the solution obtained depends critically on the initial guess of the decision variables for the optimization process. Furthermore, for the calculation of
PT
the gradient, it requires to solve the dynamic system forward in time and the corresponding co-state system backward in time. Furthermore, it is required
CE
to calculate the integration of the gradient of the Hamiltonian function over the planning horizon. Although all these tasks are to be carried out numerically, these tasks are known to be rather time consuming. In many real
AC
industrial applications, only input-output data are available [15]. For these situations, it is impossible to derive the gradient formulas. Thus, it is required to seek derivative-free optimization algorithms for finding the controls and the terminal time.
5
ACCEPTED MANUSCRIPT
In this paper, we propose to adopt a numerical computational approach based on particle swarm optimization (PSO) algorithm. PSO [16] is a derivative-
CR IP T
free and population-based stochastic search algorithm, which is easy to implement and has been found to be robust and fast in solving many practical
problems in sciences and engineering [17]. It is known that these PSO-based optimization algorithms can get trapped in local minima during the later stage of the searching process. Therefore, extensive research has been car-
AN US
ried out by many researchers over the years for overcoming this drawback. Improvements achieved are basically in the following three aspects: topology design [18], parameter adjustment [19], and hybrid algorithm with other optimization algorithms [20]. Many improved particle swarm optimization algorithms can enhance the performance of the original method at the expense
M
of increasing the complexity of these algorithms. For the topology design, it determines the social information sharing mechanism among populations.
ED
Subsequently, it changes the effectiveness of PSO. With this motivation, in this paper, a basic PSO with fully informed topology proposed in [21] is uti-
PT
lized. It is well recognized that this topology is simple and is more robust than other topologies.
CE
The rest of the paper is organised as follows. Problem formulation is given in Section 2. The problem approximation through the application of the control parameterization technique is detailed in Section 3. In Section 4,
AC
an effective optimization method based on a fully-informed particle swarm optimization method is presented and two numerical examples are solved. Finally, some concluding remarks are made in Section 5.
6
ACCEPTED MANUSCRIPT
2. Problem statement Consider the following time delayed system t ∈ [0, T ],
t ≤ 0,
x(t) = φ(t),
(1)
CR IP T
˙ x(t) = f (t, x(t), x(t − α), u(t)),
(2)
where T > 0 is an unknown terminal time; x(t) = [x1 (t), . . . , xn (t)]> ∈ Rn
is the state vector; x(t − α) ∈ Rn is the delayed state vector; α is a given
AN US
time delay; u(t) = [u1 (t), . . . , um (t)]> ∈ Rm is an unknown control function. f = [f1 , . . . , fn ]> and φ(t) = [φ1 , . . . , φn ]> are given functions.
A control function u = [u1 , . . . , um ]> is said to be an admissible control if the following boundedness constraints are satisfied i = 1, . . . , m,
t ∈ [0, T ],
(3)
M
bi ≤ ui (t) ≤ ci ,
where 0 ≤ bi < ci , i = 1, . . . , m, and bi , ci , i = 1, . . . , m, are, respectively,
ED
the lower and upper bounds of the ith control variable ui (t). Let U be the class of all admissible controls.
PT
Furthermore, we impose the following constraint: Tmin ≤ T ≤ Tmax
(4)
CE
where Tmin and Tmax are the lower and upper bounds of the terminal time, respectively.
AC
We assume that the following conditions are satisfied.
Assumption 1. The function f is continuously differentiable with respect to x and u for each t ∈ [0, T ], and piecewise continuous with respect to t for
each (x, u) ∈ Rn × Rm . Furthermore, the function φ is twice continuously differentiable. 7
ACCEPTED MANUSCRIPT
Assumption 2. There exists a real number L1 > 0 such that ˆ , u)| ≤ L1 (1 + |x| + |ˆ |f (t, x, x x| + |u|),
ˆ , u) ∈ R × Rn × Rn × Rm , (t, x, x
CR IP T
where | · | denotes the Euclidean norm.
On the basis of Assumptions 1 and 2, the dynamic system (1)-(2) admits a unique solution corresponding to each control u satisfying (3)[22]. Let
u ∈ U.
AN US
x(·|u) denote the solution of the dynamic system (1)-(2) corresponding to
We further require that the following state constraints are satisfied. Φl (t, x(t|u)) ≥ 0,
l = 1, . . . , r,
t ∈ [0, T ],
(5)
where the functions Φl , l = 1, . . . , r, are continuously differentiable with
M
respect to x for each t ∈ [0, T ] and piecewise continuous with respect to t for each x ∈ Rn . as follows:
ED
A free terminal time optimal control problem can now be formally stated
PT
Problem (P). Given system (1)-(2), find a terminal time T and an admis-
CE
sible control u(t) such that the performance function Z T J(T, u) = Φ0 (T, x(T |u)) + L0 (t, x(t|u), u)dt,
(6)
0
AC
is minimized subject to constraints (3)-(5). 3. Problem approximation In this section, we shall apply the control parameterization technique [23] to approximate Problem (P) by a sequences of nonlinear programming 8
ACCEPTED MANUSCRIPT
problems. The control parameterization technique approximates the control by a piecewise-constant function or piecewise smooth function with possible
CR IP T
discontinuities at the N − 1 pre-assigned partition time points, where the
heights of the piecewise-constant function are regarded as decision variables. Let uN i (t)
=
N X
σiN,k χ[τk−1 ,τk ) (t),
t ∈ [τk−1 , τk ),
k=1
AN US
Here, χ[τk−1 ,τk ) (t) is the indicator function defined by 1, if t ∈ [τk−1 , τk ), χ[τk−1 ,τk ) (t) = 0, elsewhere
i = 1, . . . , m.
(7)
(8)
and τk , k = 1, . . . , N , are the control switching times satisfying 0 = τ0 < τ1 < · · · < τN −1 < τN = T . For each k = 1, . . . , N and i = 1, . . . , m,
M
σiN,k is the control height of the ith piecewise constant control on the time
interval [τk−1 , τk ). Denote σ N = [(σ N,1 )> , . . . , (σ N,N )> ]> , where σ N,k =
ED
N,k > ] , k = 1, . . . , N . [σ1N,k , . . . , σm
Let ∆τk = τk − τk−1 , where the convention of ∆τ0 = 0 is adopted. It is
CE
PT
assumed that mink=1,...,N ∆τk > 0. Clearly, τk =
k X
∆τk .
i=1
In the classical control parameterization method, ∆τk is fixed. In this
AC
case, only the control heights are decision variables. Thus, to obtain an accurate optimal control policy, the number of N is required to be very large, resulting in a large number of decision variables. Thus, a non-uniform control vector parameterization is introduced in [24, 25], where the time durations ∆τk ,k = 1, ..., N , are also taken as decision variables. Let ∆τ N = 9
ACCEPTED MANUSCRIPT
[∆τ1 , . . . , ∆τN ]> . Then ∆τk and σiN,k , k = 1, . . . , N , i = 1, . . . , m, are decision variables. Let T N denote the set of all such ∆τ N and let ΞN be the set
CR IP T
containing all such σ N . Substituting (7) into system (1)-(2) gives the following system N X k=1
f (t, x(t), x(t − α), σ N,k )χIk (t),
x(t) = φ(t),
t ∈ Ik ,
t ≤ 0,
k = 1, . . . , N
(9)
(10)
P Pk−1 ∆τi , ki=1 ∆τi ). where Ik = [ i=1
AN US
˙ x(t) =
Note that, for each pair (σ N , ∆τ N ) ∈ ΞN × T N , the dynamic system
(9)-(10) admits a unique solution, which is denoted by x(·|σ N , ∆τ N ). For the constraints (4)-(5), they become
k=1
∆τk ≤ Tmax ,
M
Tmin ≤
N X
ED
bi ≤ σik ≤ ci , i = 1, . . . , m,
Φl (t, x(t|σ N , ∆τ N ) ≥ 0,
k = 1, . . . , N,
l = 1, . . . , r,
t ∈ [0, τN ].
(11) (12) (13)
PT
The new objective function is
CE
˜ N , ∆τ N ) = Φ0 (τN , x(τN |σ N , ∆τ N )) J(σ N Z τk X L0 (t, x(t|σ N , ∆τ N ), σ N,k )dt. + k=1
(14)
τk−1
AC
We may now state the approximate optimal control problem with contin-
uous state inequality constraints as follows: Problem (PN ). Given system (9)-(10), find a time duration vector ∆τ N
and a control parameter vector σ N such that the objective function (14) is minimized subject to constraints (11)-(13). 10
ACCEPTED MANUSCRIPT
Note that unless the dynamic system (9)-(10) is linear and the objective function (14) is linear or quadratic, they are non-convex. Thus, Problem (PN )
CR IP T
is non-convex. For example, we consider the optimal harvesting problem in [27]. Let us just discuss the situation when the partition number N = 1. The ˜ 1 , ∆τ1 )) three-dimensional drawing with respect to each pairs of (σ 1 , ∆τ1 , J(σ are shown in Figure 1, where the pairs of (σ 1 , ∆τ1 ) corresponding to the white
blocks are not feasible solutions. Thus, the feasible domain of Problem (PN )
AN US
is non-convex. In the following section, we will adopt an optimization strat-
egy developed based on a particle swarm optimization algorithm to solve Problem (PN ).
−5
M
−10 −15 −20
ED
The value of the cost function J
0
−25 1
Figure
1:
PT
u(t) 2
The
3
three-dimensional
4
6
drawing
t
2
with
respect
0
to
each
pairs
of
˜ 1 , ∆τ1 )) (σ , ∆τ1 , J(σ
AC
CE
1
4. Problem solving 4.1. A fully-informed particle swarm optimization Particle swarm optimization (PSO) is inspired by the social and cooperative behaviour of flocks of birds searching for food. A potential solution 11
ACCEPTED MANUSCRIPT
is called particle. A population of potential solutions is called swarm. The swarm is initialised with randomly generated particles, and the search for an
CR IP T
optimal result is carried out by updating the particles’ velocities and positions at each iteration. Recently, more and more efforts are devoted to the
design of methods for updating velocities and positions so as to avoid getting trapped in a local optimum or losing important information. In [21],
a fully informed particle swarm optimization (FIPSO) is proposed. In a
FIPSO are updated as follows. → − − → → − − v j+1 = χ(→ v ji + ϕ( P m − X ji )), i − → − → − = X ji + → v j+1 , X j+1 i i
AN US
D-dimensional search space, the velocities and positions of the particles in
i = 1, . . . , M,
j = 1, . . . , M axiter (15) (16)
i
M
where M is the population size; M axiter is the maximum iteration number; − → → − v j ∈ RD and X j ∈ RD are the velocity and position of the ith particle in the i
ED
population at the jth iteration, respectively; χ = 0.7298 is the constriction → − coefficient for the velocity; ϕ = 4.1 is the acceleration coefficient; P m refers
AC
CE
PT
to the comprehensive information learned from the ith particle’s neighbours. → − In this deterministic model, P m is calculated as → − Pm =
P
N→ − − W(k)→ ϕk Pk P , → − k∈N W(k) ϕ k → − → − ϕ k = U 0, ϕmax /|N | , k∈N
(17) (18)
→ − where N is the best among the neighbors of the particle; P k is the best position found by individual k, k ∈ N ; W(k) is the fitness difference between
the best position of the particle and the current individual, or it is a given
12
ACCEPTED MANUSCRIPT
constant. Furthermore, the velocity and position must satisfy i = 1, . . . , M,
j = 1, . . . , M axiter
(19)
i = 1, . . . , M,
j = 1, . . . , M axiter
(20)
CR IP T
− v min ≤ → v ji ≤ v max , − → X min ≤ X ji ≤ X max ,
where v min and v max are, respectively, the lower and upper bounds of the
velocity; X min and X max are, respectively, the lower and upper bounds of the position.
AN US
In the fully informed neighborhood, all neighbors are considered as sources of influence. Thus, the size and the topology of neighbors determine how diverse the influence will be. Experimental results in [21] show that FIPSO with the square topology, where the reference to the particle’s own index is removed from the neighborhood, converges faster than FIPSO with other
M
topologies. Furthermore, it has better chance of finding the best optimum.
ED
Thus, we will use FIPSO to solve the approximate optimal control problem. 4.2. A computational method based on FIPSO
PT
It is known that PSO can only solve problems with bound constraints, such as constraints (11) and (12). For optimization problems with complex
CE
constraints, such as continuous state inequality constraints, they cannot be
AC
solved directly by PSO. Thus, we introduce the following penalty function
13
ACCEPTED MANUSCRIPT
[26] to handle these complex constraints.
+
r X l=1
ρl
Z
τN
N Z X k=1
τk−1
L0 (t, x(t|σ N , ∆τ N ), σ N,k )dt
min{Φl (t, x(t|σ N , ∆τ N )), 0}2 dt
0
τk
CR IP T
J˜ρ (τN , σ N , ∆τ N ) = Φ0 (τN , x(τN |σ N , ∆τ N )) +
+ ρr+1 min{τN − Tmax , 0}2 + min{Tmin − τN , 0}2
(21)
AN US
where ρl , l = 1, . . . , r + 1, are given large penalty factors. The last three terms in (21) represent the constraint violations. Problem (PN ) is now approximated as:
Problem (PN,ρ ). Given system (9)-(10), find a time duration vector ∆τ N and a control vector σ N such that the objective function (21) is minimized
M
subject to bound constraint (12).
ED
Clearly, Problem (PN ) is a more accurate approximation of Problem (P ) when the partition number N increases. By using arguments similar to those given in [23], it can be shown that the approximate optimal cost will converge
PT
to the true optimal cost as N → ∞. However, in practice, we cannot let N = ∞. Furthermore, when N increases the decision parameters of the
CE
approximated problem will grow by m times. From our extensive numerical experiments, it suggests that the improvement in the cost value for N larger
AC
than N = 20 appears to be highly insignificant. Thus, we start with a small N , then Problem (PN,ρ ) is solved by FIPSO to obtain the optimal control.
We then increase the value of N and re-calculate the optimal control of the corresponding Problem (PN,ρ ). We repeat this process until the reduction in
the cost value is negligible. Details are listed as follows. 14
ACCEPTED MANUSCRIPT
Step 1: Initialize N . Set the swarm size M , the neighbour particles for each particle, and the maximum iteration or terminal conditions of FIPSO.
CR IP T
Step 2: Initialize the velocity and position of each particle in the swarm and its neighbours. Let the iteration p = 1. Denote the swarm by Swarmp .
Step 3: For each particle in the swarm, execute the following step: Solve the state differential equation (9)-(10) with the initial condition forward
AN US
in time from t = 0 to t = τN . Let the solution obtained be denoted as x(·|σ N , ∆τ N ). Substitute x(·|σ N , ∆τ N ) into the objective function (21) to obtain the cost value for each particle.
Step 4: Update each particle’s personal best position and the global best of the swarm. If the terminal conditions of FIPSO are satisfied, then go to
M
Step 6, else let p := p + 1 and go to Step 5.
ED
Step 5: Update the velocities and positions of the swarm according to (15)-(16) to obtain a new swarm Swarmp+1 . Go to Step 3.
PT
Step 6: Obtain the best control parameter vector σ N ∗ , time duration vector ∆τ N ∗ and the minimum cost value J˜ρ (τ N,∗ , σ N,∗ , ∆τ N,∗ ), and the best
CE
swarm P ∗M . If J˜ρ (τ N,∗ , σ N,∗ , ∆τ N,∗ ) − J˜ρ (τ N +1,∗ , σ N +1,∗ , ∆τ N +1,∗ ) ≤ ,
where = 10−6 is a required precision, then stop. Else let N := N + s,
AC
where s is some positive integer. Go to Step 2.
15
ACCEPTED MANUSCRIPT
5. Numerical results 5.1. Example 1
CR IP T
The first example involves optimal fishery harvesting, where the growth
of the population has limited resources (food, water, space). Suppose the growth of the population is limited by a carrying capacity K. Let x be the
current fish population and r be the intrinsic growth rate. The dynamics
of the population from generation to generation is governed by a logistic x(t) ]. K
Moreover, considering that the
AN US
model in the form of x(t) ˙ = rx(t)[1 −
juvenile fish cannot be harvested, and they needs time to reach maturity for reproduction, the population growth model should contain time delay from their birth to maturity. Hence, the following delayed logistic model with
M
harvesting for fish growth is considered [27]:
x(t) ˙ = 3x(t)[1 − 0.2x(t − 0.5)] − u(t), t ≤ 0,
(22) (23)
ED
x(t) = 2,
t ∈ [0, T ]
where u(t) is the harvesting effort, r = 3, K = 5.
PT
Clearly, the fish population increases, initially. However, when it gets closer to the carrying capacity, it becomes saturated forcing the rate of growth
CE
to decrease. To maximize sustainable yield and to achieve the commercial purpose of the fishery, it is necessary to harvest the population at appropriate
AC
time points with suitable harvest effort in a long run. Hence, the following control objective function are proposed : Z T Φ0 (T, u(t)) = 0.1T 2 + e−βt (CE x(t)−1 u(t)3 − pu(t))dt 0
16
(24)
ACCEPTED MANUSCRIPT
where β = 0.05 is the discount rate of the economic benefit influenced by the market; CE = 0.2 is the harvesting cost, p = 2 is the fish price. Over all, the
CR IP T
first term of equation (24) is used to shorten the total investment time. The second term is expenditure minus income. To minimize the second term is to maximum the economic benefit. In addition, to avoid overfishing and avoid exceeding the carrying capacity limitation, the following constraints must be satisfied: t ∈ [0, T ]
AN US
2 ≤x(t) ≤ 4,
(25)
0 ≤u(t) ≤ 3.5
(26)
0.1 ≤ T ≤ 100
(27)
The optimal fishery harvesting problem may now be described as follows.
M
Given the time delayed system (22)-(23), find a terminal time T and a harvesting strategy u such that the objective functional (24) is minimized,
ED
subject to constraints (25)-(27). By the approximation method described in Section 3, the free terminal
PT
time optimal control problem is approximated by the following optimization problem. Given the dynamic system
CE
N X
˙ x(t) =
AC
k=1
x(t) = 2,
{3x(t)[1 − 0.2x(t − 0.5)] − σ k }χIk ,
t ∈ Ik ,
k = 1, . . . , N (28)
t ≤ 0,
(29)
P Pk N where Ik = [ k−1 and a control i=1 ∆τi , i=1 ∆τi ), find a positive vector ∆τ 17
ACCEPTED MANUSCRIPT
parameter vector σ N such that the cost functional N
N Z X k=1
τk
τk−1
e−0.05t [0.2x(t)−1 (σ k )3 − 2σ k ]dt
CR IP T
J˜ρ (τN , σ , ∆τ ) = N
N N nX o2 X + 0.1 ∆τk + ρ1 (min{ ∆τk − 0.1, 0})2 k=1
k=1
+ ρ1 (min{100 − τN
0
k=1
∆τk , 0})2
[(min{x(t) − 2, 0})2 + (min{4 − x(t), 0})2 ]dt
AN US
+ ρ2
Z
N X
(30)
is minimized subject to the bound constraints (11) and (12), where ρ1 and ρ2 are penalty factors.
M
We solve this problem by writing a Matlab program on a personal computer with the following configuration: Intel Core i5-4590 3.3GHz CPU, 8G
ED
RAM, 64-bit Window 7 operating system. In order to test the performance of the proposed method, we choose N = 9 as given in [27]. The penalty factors are ρ1 = 1000 and ρ2 = 105 . The performances of the FIPSO are evaluated
PT
by comparing with that of the comprehensive learning particle swarm optimizer (CLPSO) proposed in [28], genetic algorithm (GA) given in [29]. For
CE
all three optimization algorithms, the swarm size or population size are 30, and the maximum iteration is 1000. According to reference [21], for FIPSO,
AC
the square topology is used for choosing four neighbours, and W(k) is the fitness difference between the best position of the particle and the current individual. For GA, the selection rate is 0.9, cross rate is 0.8, and mutation rate is 0.1. The other parameters of CLPSO are as defined in [28]. The optimal results are summarized in Table 1 and Figure 2-Figure 4. 18
ACCEPTED MANUSCRIPT
Table 1: Comparisons results of Example 1
FIPSO
CLPSO
GA
reference [14]
optimal cost value
-26.032
-26.025
-25.620
-25.97
optimal terminal time
12.1221
12.1487
12.1448
12.1565
CR IP T
optimization algorithm
AN US
Table 1 shows that FIPSO solves the optimal control problem with a shorter
terminal time (12.12 years), and a smaller cost value. Compared with [27], the harvest time can be shorted by 10 days. Figure 2 demonstrates that the convergence speed of FIPSO is faster than that of CLPSO and GA optimization algorithm. And the optimal state and control of FIPSO have a tendency
M
of becoming stable. Hence, FIPSO is superior to the other two optimization
are satisfied. 10
x 10
4
FIPSO
cost function value
CE AC
CLPSO
GA
−24.5
PT
8
ED
methods in solving this optimal control problem. In addition, all constraints
6
−25
4
−25.5 −26 490
2
495
500
0 100
200
300
400 500 600 iteration number
700
800
900
1000
Figure 2: The convergence trajectories of FIPSO, GA, and CLPSO.
19
ACCEPTED MANUSCRIPT
CR IP T
4
x(t)
3.5 3
CLPSO FIPSO GA
2.5
0
2
4
6
8
t
10
12
AN US
2
M
Figure 3: The optimal state trajectories obtained by FIPSO, GA, and CLPSO, respectively
ED
3.5
PT
u(t)
3
2.5
FIPSO CLPSO GA
1.5
0
2
4
6
AC
CE
2
t
8
10
12
Figure 4: The optimal piecewise-constant controls obtained by FIPSO, GA, and CLPSO, respectively
20
ACCEPTED MANUSCRIPT
Table 2: Comparisons results of Example 2
CLPSO
GA
reference [1]
CR IP T
optimization algorithm FIPSO optimal cost value
5.287
5.425
8.875
5.37
optimal terminal time
1.349
1.1818
1.005
1.446
AN US
5.2. Example 2
The second example is the optimization problem given in [1]
min J = subject to the dynamics
T
0
(x2 (t) + u2 (t))dt + 0.1T
M
u,T
Z
ED
x(t) ˙ = x(t) + 0.8x(t − 0.15) + 0.9u(t),
t ∈ [0, T ]
(31)
(32)
t ≤ 0,
(33)
−10 ≤ u(t) ≤ 10,
(34)
PT
x(t) = 1,
We solve this problem with control partition numbers N = 6, using
CE
FIPSO, CLPSO, GA algorithms with the same parameter settings for solving example 1. Comparison results are given in Table 2. We can also see that
AC
FIPSO is more effectiveness than the other optimization algorithms. From the results given in Table 2, we can see that the minimum cost has been reduced by 1.6% and the minimum time has been reduced by 7% compared with that of [1].
21
ACCEPTED MANUSCRIPT
6. Conclusion In this paper, a class of free terminal time optimal control problem in-
CR IP T
volving time delayed system is investigated. To solve this problem, we apply the control parameterization method to yield an approximate optimization
problem. Then a simple effective optimization method is developed, based on a fully informed particle swarm optimization method, to solve the approx-
imate problem. Numerical results show that the proposed method is more
AN US
effective than other existing solving methods with shorter terminal time and better cost function values. Acknowledgment
M
The authors thank the anonymous reviewers and Mr. Jianhua Liu for their valuable comments and suggestions to improve this paper. This work
ED
was supported by the Natural Science Foundation of Fujian Province, China (Grant No. 2016J05154), the National Natural Science Foundation of China
PT
(Grant No. 61673116, 61503131).
CE
References
[1] A. Boccia, P. Falugi, H. Maurer, R. B. Vinter. Free time optimal control problems with time delays. IEEE Conference on Decision and Control,
AC
2013:520-525.
[2] M. Liu, C. Bai. Optimal harvesting of a stochastic logistic model with time delay. Journal of Nonlinear Science, 2015, 25(2):277-289.
22
ACCEPTED MANUSCRIPT
[3] C. Liu, Z. Gong. Modelling and optimal control of a time-delayed switched system in fed-batch process. Journal of the Franklin Institute,
CR IP T
2014, 351(2): 840-856. [4] I. Elmouki. Free terminal time optimal control problem for the treatment of HIV infection. Jawra Journal of the American Water Resources Association, 2016, 7(1):148-161.
AN US
[5] J. J. Ewoud, Q. Chu, G. C. D. Croon. Adaptive incremental nonlinear
dynamic inversion for attitude control of micro air vehicles.Journal of Guidance Control and Dynamics, 2016, 39(3): 450-461. [6] L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, E. F.
science, 1962.
M
Mishchenko, The mathematical theory of time optimal processes, Inter-
ED
[7] R. K. Biswas, S. Sen. Free final time fractional optimal control problems. Journal of the Franklin Institute, 2014, 351(2):941-951.
PT
[8] S. Pooseh, R. Almeida, D. F. M. Torres. Fractional Order Optimal Control Problems with Free Terminal Time. Journal of Industrial and Man-
CE
agement Optimization, 2014, 10(2): 363-381. [9] J. Yuan, C. Liu, X. Zhang, J. Xie, E. Feng, H. Yin, Z. Xiu. Optimal
AC
control of a batch fermentation process with nonlinear time-delay and free terminal time and cost sensitivity constraint. Journal of Process Control, 2016, 44: 41-52.
[10] R. Li, Y. Shi, Finite-time optimal consensus control for second-order 23
ACCEPTED MANUSCRIPT
multi-agent systems. Journal of Industrial & Management Optimization, 2014, 10(3):929-943.
CR IP T
[11] C. Jiang, Q. Lin, C. Yu, K. L. Teo, G. R. Duan. An exact penalty method
for free terminal time optimal control problem with continuous inequality constraints. Journal of Optimization Theory and Applications, 2012, 154(1):30-53.
AN US
[12] Z. Gong, C. Liu, Y. Wang. Optimal control of switched systems with multiple time-delays and a cost on changing control. Journal of Industrial and Management Optimization, 2017, 13, doi:10.3934/jimo.2017042.
[13] J. M. Jeong, S. J. Son, Time optimal control of semilinear control sys-
M
tems involving time delays. Journal of Optimization Theory and Appli-
ED
cations, 2015, 165(3):793-811.
[14] C. V. Peroni, N. S. Kaisare, J. H. Lee. Optimal control of a fed-batch bioreactor using simulation-based approximate dynamic programming.
CE
790.
PT
IEEE Transactions on Control Systems Technology, 2005, 13(5): 786-
[15] T. J. J. V. D. Boom, J. B. Klaassens, R. Meiland. Real-time time-
AC
optimal control for a nonlinear container crane using a neural network. 2007:79-84.
[16] R. Eberhat, J. Kennedy, A new optimizer using particle swarm theory. Proceedings of 6th International Symposium On Mico Machine and Human Science, 1995: 39-43. 24
ACCEPTED MANUSCRIPT
[17] C. Ou, W. Lin. Comparison between PSO and GA for parameters optimization of PID controller. IEEE International Conference on Mecha-
CR IP T
tronics and Automation. 2006: 2471-2475. [18] Q. Liu, W. Wei, H. Yuan, Z. H. Zhan, Y. Li. Topology selection for particle swarm optimization. Information Sciences,2016, 363: 154-173.
[19] A. B. Hashemi, M. R. Meybodi. A note on the learning automata based
puting, 2011, 11(1):689-705.
AN US
algorithms for adaptive parameter selection in PSO. Applied Soft Com-
[20] H. Garg. A hybrid PSO-GA algorithm for constrained optimization problems. Applied Mathematics & Computation, 2016, 274(11):292-305.
M
[21] R. Mendes, J. Kennedy, O. Neves, The fully informed particle swarm: simpler, maybe better. IEEE Transactions on evolutionary computation,
ED
2004, 8(3):204-210
[22] N. U. Ahmed, Dynamic systems and control with applications, Singa-
PT
pore: World Scientific, 2006. [23] K. L. Teo, C. J. Goh, K. H. Wong, A unified computational approach
CE
to optimal control problems, Essex: Longman Scientific and Technical,
AC
1991.
[24] Y. Lei, S. Li, Q. Zhang, X. Zhang, Free final time optimal control based on non-uniform control vector parameterization. Journal of Systems Science and Mathematical Sciences, 2012, 32(3): 277-287.
25
ACCEPTED MANUSCRIPT
[25] Q. Lin , R. Loxton, K. L. Teo, The control parameterization method for nonlinear optimal control: A survey. Journal of Industrial & Manage-
CR IP T
ment Optimization, 2017, 10(1):275-309. [26] C. J. Yu, K. L. Teo, L. S. Zhang, Y. Q. Bai, On a refinement of the con-
vergence analysis for the new exact penalty function method for continuous inequality constrained optimization problem. Journal of Industrial
AN US
Management and Opti- mization, 2012, 8(2):485-491.
[27] C. Liu, R. Loxton, K. L. Teo. A computational method for solving timedelay optimal control problems with free terminal time. Systems & Control Letters, 2014, 72: 53-60.
[28] J. J. Liang, A. K. Qin, P. N. Suganthan, S. Baskar, Comprehensive learn-
M
ing particle swarm optimizer for global optimization of multimodal func-
281-295.
ED
tions. IEEE Transactions on Evolutionary Computation, 2006, 10(3):
PT
[29] Goldberg D E. Genetic algorithm in search, optimization, and machine
AC
CE
learning. 1989, xiii(7): 2104C2116.
26