# Time Optimal Control of a Thermoelastic System

Authors: Qiyu Liu; Qunxiong Zhu; Zhiqiang Geng; Longjin Lv
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This paper considers the numerical approximation for the time optimal control problem of a thermoelastic system with some control and state constraints. By the Galerkin finite element method (FEM), the original problem is projected into a semidiscrete optimal control problem governed by a system of ordinary differential equations. Then the optimal time and control parameterization method is applied to reduce the original system to an optimal parameter selection problem, in which both the optimal time and control are taken as decision variables to be optimized. This problem can be solved as a nonlinear optimization problem by a hybrid algorithm consisting of chaotic particle swarm optimization (CPSO) and sequential quadratic programming (SQP) algorithm. The numerical simulations demonstrate the effectiveness of the proposed numerical approximation method.


---

## Body

## 1. Introduction

The expansion and contraction of a material exposed to temperature changes are of great importance in many applications such as stabilization of satellite antennas, which is commonly with large temperature variations and characterized by the thermoelastic coupling [1]. This phenomenon is modelled by a system consisting of a hyperbolic elasticity equation describing the displacement coupled with a parabolic equation for the temperature. Over the past few decades, there has been much work on control of thermoelastic systems including [1–8] and references therein. For more studies of the related problems of the thermoelasticity, we would like to mention some other related works [9, 10]. Very little, however, is known about the time optimal control of thermoelastic systems.As one of the most important optimal control problems, time optimal control problems aim at driving the corresponding trajectory of the system from an initial state to a given target set in the shortest time, through applying constrained controllers. Time optimal control problems in the finite-dimensional case [11, 12] were first studied, which were governed by a system of ordinary differential equations. Then time optimal control problems were extended to the infinite-dimensional case; namely, the controlled equations are partial differential equations (PDEs) or abstract evolution equations in Banach spaces [13–16]. The main theoretical tool to deal with time optimal control problems analytically is the Pontryagin minimum principle [17]. However, practical problems are usually becoming too complex to allow analytical techniques alone. Thus numerical methods for solving these problems attracted more and more researchers’ attention [18–22].In this paper, we will study the numerical approximation for the time optimal control problem governed by the one-dimensional linear thermoelastic system. As far as we know, little work has been done on numerical methods for solving this problem. Our method follows the discrete-then-optimize approach [23–25], whereby the PDE system is first projected onto the finite-dimensional subspace to obtain a system of ordinary differential equations (ODEs), and computational optimal control techniques are applied; then the resulting problem is solved by appropriate optimization methods. Specifically, by making use of the standard finite element method (FEM), the original problem was firstly projected into a time optimal control problem governed by a system of ODEs. Noting that the terminal time of this problem is unspecified, we will, then, introduce a transformation such that the optimal time can be regarded as an unknown system parameter. Moreover, control parametrization method, which involves approximating the control function by a piecewise-constant function with possible discontinuities at a set of preassigned switching points [26, 27], will be carried out and the original problem can be reduced to an optimal parameter selection problem. In principle, this problem can be solved as a nonlinear optimization problem by standard mathematical programming algorithms such as SQP.The SQP method is one of the nonlinear programming methods for constrained optimization. The basic idea of SQP is to make a quadratic approximation of the Lagrange function of the original problem, to form a quadratic programming (QP) subproblem. The solution of QP subproblem is then used to form a search direction for a line search procedure in the next iteration. Although with the strong capability of local searching, the SQP method is dependent on the initial estimation and is easy to be trapped into a local minimum for multimodal problem. Thus a hybrid PSO is introduced to overcome this drawback.The PSO algorithm is one of the most well-known bioinspired algorithms used in optimization problems, which basically consists of a number of particles that collectively move in search of the global optimum. Meanwhile, similar with other evolution algorithms, the traditional PSO is also prone to get trapped into the local minimum when solving some complex problems. There is a rich literature devoted to the improvement of PSO; for example, refer to [28–30] and references therein. In our work, a chaotic method is integrated with PSO to further promote the diversity of PSO [31–33]. This chaotic PSO algorithm is referred to as CPSO. In this way, the possibility of exploring a global optimum solution in problems is increased. Furthermore, as empirical and theoretical studies show that PSO is sensitive to the control parameters, especially the inertia weight and acceleration coefficients, we also make full use of existing results of parameter controlling strategies to improve the performance of CPSO.In this paper, a hybrid algorithm consisting of CPSO and SQP is introduced to solve the optimal parameter selection problem, where CPSO can be viewed as the global optimizer while the SQP is employed for the local search.The rest of this paper is organized as follows. The formulation of the time optimal control problem is introduced in Section2. In Section 3, we obtain the semidiscrete form of the original optimal control problem with FEM method. In Section 4, an optimal parameter selection problem is obtained by the optimal time and control parameterization method, and a hybrid optimization algorithm CPSO-SQP is employed to solve this problem. In Section 5, the numerical simulations are presented to illustrate the effectiveness of our numerical approximation method. Finally, in Section 6, we conclude the paper with summary comments and suggestions for future research.


## 2. Formulation of the Time Optimal Control Problem

Consider a linear one-dimensional thermoelastic rod of lengthl with unit reference density. The controlled thermoelastic system takes the form(1)wtt-awxx+γθx=f,0<x<l,0<t<∞,θt+bwxt-kθxx=g,0<x<l,0<t<∞,w0,t=θ0,t=0,wxl,t=u1t,θxl,t=u2t,0<t<∞,wx,0=w0x,wtx,0=w1x,θx,0=θ0x,where w(x,t) is the displacement, θ(x,t) represents the temperature deviation from the reference temperature of the bar at (x,t), w0, w1, and θ0 are the given initial data, f(x,·),g(x,·)∈L2(0,l) are the external effects acting on this system, and u1 and u2 are the boundary control functions. In the above, the first equation is the momentum equation and the second equation is the energy equation. We assume that the system parameters a, b, γ, and k are given constants with a,k>0, and b,γ≠0, which are dependent on the material properties.DefineV=w∈H10,l∣w0=0 with the norm (2)wV=∫0lw′x2dx1/2,∀v∈V.By the famous Poincare’s theorem, we see that V is a Hilbert space. Moreover, let H=L2(0,l) with the standard norm ∥·∥. Define H=V×H×H equipped with the norm (3)yH=∫0lwxx2dx+∫0lvx2dx+∫0lθx2dx1/2,∀y=w,v,θ⊤∈H.Then, H is a Hilbert space and the energy function of the thermoelastic system (1) is(4)Et=w,wt,θ⊤H,where (w,θ) is the solution of (1).Define(5)U=U·=u1·,u2·T∈L∞0,∞×L∞0,∞∣αi≤uit≤βia.e.t∈0,∞,i=1,2,where αi,βi∈R,i=1,2, are given constants. Let r∈R be positive constant and assume that, for given (w0,w1,θ0)∈H, there exists at least one control U∈U, such that the solution (w,θ) of (1) satisfies that E(t)≤r for some t>0.Under this setting, the time optimal control problem of the thermoelastic system, which we shall study, is formulated as follows:(6)TOCPinf⁡T∣ET;U≤r,U∈U,where the energy function E(t)=E(t;U) is defined by (4) and (w,θ) is the solution of (1) with respect to (w0,w1,θ0)∈H and U∈U. Define (7)T∗≜inf⁡T∣ET;U≤r,U∈Uas the optimal time for the problem (TOCP). Moreover, if there exists U∗∈U such that E(T∗;U)≤r, we call U∗ as an optimal control for the problem (TOCP) and the solution (w∗,θ∗) of (1) corresponding to the optimal control U∗ as the optimal state. Throughout this paper, we assume that there exists a pair (T,U) with U∈U such that E(T;U)≤r. By the standard arguments [34, 35], we can get the existence of the optimal time T∗ and the optimal control U∗.


## 3. Finite Element Approximation

In this section, we will apply the Galerkin finite element method [36, 37] to discretize the original problem (TOCP), an infinite-dimensional distributed parameter time optimal control problem, into a semidiscrete approximation problem governed by a finite-dimensional lumped parameter system.In order to construct a finite-dimensional subspaceVh of the space V, we subdivide [0,l] into N subintervals (8)Ij=xj-1,xj,j=1,2,…,N,by the grid points: (9)0=x0<x1<⋯<xN-1<xN=l,where the length of each subinterval is hj=xj-xj-1,j=1,2,…,N. Then h=max⁡hj is the mesh size of this triangulation and a measure of how fine the partition is. Note that in general the mesh points xj need not be equally spaced. For simplification, we choose a uniform spacing here.Now letVh be the set of functions vh such that vh is linear on each subinterval Ij and continuous on [0,l]: (10)Vh=vh∈C0,lvhIj∈P1Ijforj=1,2,…,N,vh0=0,where P1(Ij) is the space of polynomials of degree at most 1 defined on the subinterval Ij.We observe thatVh is a N-dimensional space and Vh⊂V. Let {ϕi;i=1,2,…,N} be the basis functions of the finite-dimensional space Vh, defined by (11)ϕixj=1,ifi=j;0,ifi≠j,for i=1,2,…,N,j=1,2,…,N, where ϕi is the continuous piecewise linear function that takes value 1 at node point xj and the value 0 at other node points.Then, the finite element approximation of (1) is the L2-projection from L2(0,l) onto Vh, defined by(12)Phv,vh=v,vh,∀v∈L20,l,vh∈Vh,where (·,·) is the inner product of L20,l. Taking the inner product of (1) with the test function vh, replacing V by the finite-dimensional subspace Vh, and applying boundary conditions, we get the following semidiscrete analogue of system (1): find wh(t)=wh(·,t),θh(t)=θh(·,t)∈Vh for t∈[0,T], such that(13)wttht,vh+awxht,vxh+γθxh,vh-au1tvhl=ft,vh,θtht,vh+bwxtht,vh+kθxht,vxh-ku2tvhl=gt,vh,wh0=Phw0,wth0=Phw1,θh0=Phθ0.Sincewh,θh∈Vh=span{ϕi}, they can be written as the linear combination of the basis functions:(14)whx,t=∑j=1NXjtϕjx,θhx,t=∑j=1NZjtϕjx,where Xj(t) and Zj(t),j=1,2,…,N, are weighting functions. Moreover, we can write(15)Phw0x=∑j=1NXj0ϕjx,Phw1x=∑j=1NXj1ϕjx,Phθ0x=∑j=1NZj0ϕjx,where Xj0(t), Xj1(t), and Zj0(t),j=1,2,…,N, are weighting functions.Substituting (14) and (15) into (13) and taking vh=ϕi,i=1,2,…,N, yield that(16)∑j=1Nϕi,ϕjX¨jt+a∑j=1Nϕi′,ϕj′Xjt+γ∑j=1Nϕi,ϕj′Zjt-aϕilu1t=ft,ϕi,∑j=1Nϕi,ϕjZ˙jt+b∑j=1Nϕi,ϕj′X˙jt+k∑j=1Nϕi′,ϕj′Zjt-kϕilu2t=gt,ϕi.Define (17)Xt=XitN×1,M=ϕi,ϕjN×N,K=ϕi′,ϕj′N×N,W=ϕi,ϕj′N×N,B=ϕilN×1,Ft=ft,ϕiN×1,Gt=gt,ϕiN×1.Thus, it follows from (16) that(18)MX¨t+aKXt+γWZt-aBu1t=Ft,MZ˙t+bWX˙t+kKZt-kBu2t=Gt.Define X˙=X1 and (19)Y=XX1Z,R=0I0-aM-1K0-γM-1W0-bM-1W-kM-1K,D=00aM-1B00kM-1B,Ht=0M-1FtM-1Gt.Thus, (13) can be rewritten as follows:(20)Y˙t=RYt+DUt+Ht,0<t≤T,Y0=Y0,where (21)Y0=X0X1Z0.Note that this system is a set of initial value problems (IVPs) and is equivalent to (13). Define (22)Eht;U=Y,ΠYR3N,where (·,·)R3N denotes the inner product in R3N and (23)Π=K000M000M.Now, by the finite element method, we get a new time optimal control problem (TOCPh) governed by a system of controlled ordinary differential equations from the original problem (TOCP) as follows: (24)TOCPhinf⁡T∣EhT;U≤r,U∈U,where Y(t):[0,T]→R3N is the solution of (20).


## 4. Time Optimal Control Computation

### 4.1. Optimal Time Parameterization

Obviously, the problem(TOCPh) is a problem with unspecified terminal time. In this case, we may regard the optimal time T as an unknown system parameter. Write (25)σ0=T,and define the transformation (26)t=σ0s,s∈0,1.Let (27)Y¯s=Yσ0s=Yt,U¯s=Uσ0s=Ut.Then, (20) is equivalent to the following equation:(28)Y¯˙s=σ0RY¯s+DU¯s+Hσ0s,0<s≤1,Y¯0=Y0.Thus, the problem (TOCPh) can be transformed to the following form: (29)TOCP¯hinfσ0∣E¯h1;U¯≤r,σ0≥0,U¯∈U,where E¯h(s;U¯)=(Y¯(s),ΠY¯(s))R3N for s∈[0,1].By this optimal time parameterization process, the problem(TOCPh) becomes a combined optimal parameter selection and optimal control problem, which has a fixed time horizon [0,1].


### 4.2. Control Parameterization

To approximate the problem(TOCP¯h), we will use the control parametrization method [27, 38]. More specifically, the admissible controls will be restricted to a suitable linear combination of basis functions, where the coefficients in the linear combination are decision variables to be selected optimally. This produces a sequence of finite-dimensional optimization problems for approximating problem (TOCP¯h). As piecewise-constant basis functions are normally applied in control parameterization, the optimal control problem (TOCP¯h) can be reduced to an optimal parameter selection problem.We subdivide the time horizon[0,1] into p subintervals [sk-1,sk),k=1,2,…,p, where p≥1 is the number of approximation subintervals, [sk-1,sk) is the kth approximation subinterval, and sk,k=0,1,…,p, are fixed knot points such that(30)0=s0<s1<s2<⋯<sp-1<sp=1.Define(31)Up=σ=σ1,σ2,…,σp∈R2×p∣σk=σ1k,σ2k∈R2,αi≤σik≤βi,k=1,2,…,p,i=1,2.The vector-valued control function U(t) can be approximated by a constant vector on each subinterval: (32)Ut≈σk,fort∈tk-1,tk.Then we restrict the control function to the piecewise-constant control of the form(33)Ut≈Upt∣σ=u1pt,u2ptT=∑k=1pσkχtk-1,tkt,t∈0,t,where σ∈Up and χ[tk-1,tk) is the indictor function of the interval [tk-1,tk) defined by (34)χtk-1,tkt=1,ift∈tk-1,tk,0,otherwise.Corresponding to optimal time parameterization, the control function is transformed as follows:(35)U¯s≈Ups∣σ=∑k=1pσkχsk-1,sks,s∈0,1,where σ∈Up and χ[sk-1,sk) is the indictor function of the interval [sk-1,sk) defined by (36)χsk-1,sks=1,ifs∈sk-1,sk,0,otherwise.By taking the approximationUp(s)=∑k=1p‍σkχ[sk-1,sk)(s) of U¯, (28) takes the following form:(37)Y¯˙s=σ0RY¯s+D∑k=1pσkχsk-1,sks+Hσ0s,0<s≤1,Y¯0=Y0.Let (38)U¯p=σ¯=σ0σ=σ0,σ1,σ2,…,σp⊤∈Rp+1∣σ0≥0,σ∈Upand Yp·∣σ¯ denote the solution of system (37) corresponding to σ¯∈U¯p.Now, we may state the approximate optimal parameter selection problem as follows:(39)TOCP¯phinf⁡σ0∣Eph1∣σ¯≤r,σ¯∈U¯p,where Eph(s;σ¯)=(Y¯p(s),ΠY¯p(s))R3N and Y¯p(s)=Y¯ps∣σ¯:[0,1]→R3N is the solution of (37).


### 4.3. Optimization Implementation

The finite element method and the control parameterization method have been applied to various optimal control problems. The convergence for approximation schemes of these two methods has been discussed by many authors (see [18, 20, 39] for the FEM or [26, 38, 40] for control parameterization method). Inspired by these results, by making use of the classical space discretization scheme with FEM, we first project the original time optimal control problem (TOCP) into a semidiscrete problem (TOCPh) governed by a system of ODEs. Then, through the time parameterization and control parameterization, the semidiscrete approximation (TOCPh) of the original problem has been reduced to an optimal parameter selection problem (TOCP¯ph) in the canonical form, where σ0≥0 is the time parameter and σ¯∈Up is the control parameter. The optimal parameter selection problem is basically nonlinear mathematical programming problem and hence many gradient-based optimization techniques can be applied to optimal parameter selection problem. Combined with the algorithms for computing the gradients of the cost and constraint functions such as the variation method [41, 42] and the costate method [27, 43], the optimal parameter selection problem (TOCP¯ph) can be solved numerically with some mathematical programming algorithms.In this paper, a hybrid optimization algorithm CPSO-SQP from [31, 32], which combines SQP and CPSO, is introduced to solve the problem (TOCP¯ph). In this algorithm, the CPSO is taken as the global main optimizer of the algorithm while the SQP is employed for the local search to improve the solution.In PSO, each solution is called a “particle” flying around in theD dimensional search space with a velocity Vi=(vi1,vi2,…,viD). During flight, each particle of a swarm adjusts its position Xi=(xi1,xi2,…,xiD) according to the best position encountered by itself (personal best position Pi=(pi1,pi2,…,piD) of its previous flights) and its neighbors (global best position Pg=(pg1,pg2,…,pgD) of the population) as the following formulas:(40)Vik+1=w·Vik+c1·rand·Pi-Xik+c2·rand·Pg-Xik,Xik+1=Xik+Vik+1,where i=1,2,…,M, M is the swarm size, w is the inertia weight, c1 and c2 are positive constants called acceleration coefficients for regulating the maximum step size for the particle i to fly towards Pi and Pg, respectively, superscript k represents the algorithm iteration count, and rand() is a random value in the range [0,1]. Xi and Vi should be clamped to the range of [Xmin,Xmax] and [Vmin,Vmax] to control excessive roaming of particles outside the search space, respectively.Because of its ergodicity and stochasticity, the piecewise linear chaotic map is used to improve the global search capability of PSO, which is defined as follows:(41)cxik+1=cxikq,cxik∈0,q,1-cxik1-q,cxik∈q,1,where cxik+1 denotes the chaotic mapping of ith chaotic variable in the kth iteration and q∈(0,0.5)∪(0.5,1) is a given parameter. Under this chaotic mapping, cx behaves chaotically in (0,1). Moreover, based on the fact that chaotic search is efficient in small range, we take chaotic mapping within a dynamically shrinking range around the global best position of the swarm by the following formula:(42)x^idk+1=pgd·1+η·cxik+1,where d=1,2,…,D, x^idk+1 denotes the chaotic mapping of the particle i in the dth dimension, pgd is the global best position in the dth dimension, and η is a scaling parameter with an empirical value 1.1. After chaos, the particle i should still be clamped to the range of [Xmin,Xmax] by the following formula:(43)Xik+1=max⁡min⁡X^ik+1,Xmax,Xmin.Integrated with chaotic map and SQP, the hybrid optimization algorithm CPSO-SQP is described as follows.Step 1.
Initialize swarm size, maximum function evaluations, other parameters, and the initial positionXi and velocity Vi of each particle in the swarm.Step 2.
Evaluate the objective function of each particle, then set the initial personal best solutionPi, select the global best particle Pg of CPSO to guide the swarm evolution direction, and set Po=Pg which is the global best particle of CPSO-SQP.Step 3.
Update the velocity and position of each particle according to (40) if the particle is not in stagnation, otherwise, according to (41)–(43).Step 4.
Evaluate the objective function of each particle in the new position.Step 5.
If the current position is better than personal best solution, updatePi; otherwise, set this particle in stagnation.Step 6.
IfPi is better than global best solution of CPSO, update Pg and set a flag for arousing SQP. Moreover, if Pg is better than Po, update Po.Step 7.
IfPg is just updated (the flag for SQP is true), set Pg as initial value and run SQP. If the result of SQP is better than Po, update Po.Step 8.
If maximum function evaluation is reached, outputPo as the global best solution; otherwise, go back to Step 3 for continuing iterations.In this algorithm, the powerful global search capability of CPSO to search the solution space and the solution obtained from CPSO is taken as an initial estimate of SQP algorithm to compensate for the weakness of SQP weak global search. With good local convergence performance and strong nonlinear convergence speed, SQP algorithm makes a careful search for the optimization problem to make up shortcomings of the CPSO weak local search. Thus, the possibility of exploring a global minimum in problems with more local optima is increased. Benefiting from the fast globally converging characteristics of CPSO and the effective local search ability of SQP, the proposed method can obtain the global optimal result of our problem through constant iteration.


## 4.1. Optimal Time Parameterization

Obviously, the problem(TOCPh) is a problem with unspecified terminal time. In this case, we may regard the optimal time T as an unknown system parameter. Write (25)σ0=T,and define the transformation (26)t=σ0s,s∈0,1.Let (27)Y¯s=Yσ0s=Yt,U¯s=Uσ0s=Ut.Then, (20) is equivalent to the following equation:(28)Y¯˙s=σ0RY¯s+DU¯s+Hσ0s,0<s≤1,Y¯0=Y0.Thus, the problem (TOCPh) can be transformed to the following form: (29)TOCP¯hinfσ0∣E¯h1;U¯≤r,σ0≥0,U¯∈U,where E¯h(s;U¯)=(Y¯(s),ΠY¯(s))R3N for s∈[0,1].By this optimal time parameterization process, the problem(TOCPh) becomes a combined optimal parameter selection and optimal control problem, which has a fixed time horizon [0,1].


## 4.2. Control Parameterization

To approximate the problem(TOCP¯h), we will use the control parametrization method [27, 38]. More specifically, the admissible controls will be restricted to a suitable linear combination of basis functions, where the coefficients in the linear combination are decision variables to be selected optimally. This produces a sequence of finite-dimensional optimization problems for approximating problem (TOCP¯h). As piecewise-constant basis functions are normally applied in control parameterization, the optimal control problem (TOCP¯h) can be reduced to an optimal parameter selection problem.We subdivide the time horizon[0,1] into p subintervals [sk-1,sk),k=1,2,…,p, where p≥1 is the number of approximation subintervals, [sk-1,sk) is the kth approximation subinterval, and sk,k=0,1,…,p, are fixed knot points such that(30)0=s0<s1<s2<⋯<sp-1<sp=1.Define(31)Up=σ=σ1,σ2,…,σp∈R2×p∣σk=σ1k,σ2k∈R2,αi≤σik≤βi,k=1,2,…,p,i=1,2.The vector-valued control function U(t) can be approximated by a constant vector on each subinterval: (32)Ut≈σk,fort∈tk-1,tk.Then we restrict the control function to the piecewise-constant control of the form(33)Ut≈Upt∣σ=u1pt,u2ptT=∑k=1pσkχtk-1,tkt,t∈0,t,where σ∈Up and χ[tk-1,tk) is the indictor function of the interval [tk-1,tk) defined by (34)χtk-1,tkt=1,ift∈tk-1,tk,0,otherwise.Corresponding to optimal time parameterization, the control function is transformed as follows:(35)U¯s≈Ups∣σ=∑k=1pσkχsk-1,sks,s∈0,1,where σ∈Up and χ[sk-1,sk) is the indictor function of the interval [sk-1,sk) defined by (36)χsk-1,sks=1,ifs∈sk-1,sk,0,otherwise.By taking the approximationUp(s)=∑k=1p‍σkχ[sk-1,sk)(s) of U¯, (28) takes the following form:(37)Y¯˙s=σ0RY¯s+D∑k=1pσkχsk-1,sks+Hσ0s,0<s≤1,Y¯0=Y0.Let (38)U¯p=σ¯=σ0σ=σ0,σ1,σ2,…,σp⊤∈Rp+1∣σ0≥0,σ∈Upand Yp·∣σ¯ denote the solution of system (37) corresponding to σ¯∈U¯p.Now, we may state the approximate optimal parameter selection problem as follows:(39)TOCP¯phinf⁡σ0∣Eph1∣σ¯≤r,σ¯∈U¯p,where Eph(s;σ¯)=(Y¯p(s),ΠY¯p(s))R3N and Y¯p(s)=Y¯ps∣σ¯:[0,1]→R3N is the solution of (37).


## 4.3. Optimization Implementation

The finite element method and the control parameterization method have been applied to various optimal control problems. The convergence for approximation schemes of these two methods has been discussed by many authors (see [18, 20, 39] for the FEM or [26, 38, 40] for control parameterization method). Inspired by these results, by making use of the classical space discretization scheme with FEM, we first project the original time optimal control problem (TOCP) into a semidiscrete problem (TOCPh) governed by a system of ODEs. Then, through the time parameterization and control parameterization, the semidiscrete approximation (TOCPh) of the original problem has been reduced to an optimal parameter selection problem (TOCP¯ph) in the canonical form, where σ0≥0 is the time parameter and σ¯∈Up is the control parameter. The optimal parameter selection problem is basically nonlinear mathematical programming problem and hence many gradient-based optimization techniques can be applied to optimal parameter selection problem. Combined with the algorithms for computing the gradients of the cost and constraint functions such as the variation method [41, 42] and the costate method [27, 43], the optimal parameter selection problem (TOCP¯ph) can be solved numerically with some mathematical programming algorithms.In this paper, a hybrid optimization algorithm CPSO-SQP from [31, 32], which combines SQP and CPSO, is introduced to solve the problem (TOCP¯ph). In this algorithm, the CPSO is taken as the global main optimizer of the algorithm while the SQP is employed for the local search to improve the solution.In PSO, each solution is called a “particle” flying around in theD dimensional search space with a velocity Vi=(vi1,vi2,…,viD). During flight, each particle of a swarm adjusts its position Xi=(xi1,xi2,…,xiD) according to the best position encountered by itself (personal best position Pi=(pi1,pi2,…,piD) of its previous flights) and its neighbors (global best position Pg=(pg1,pg2,…,pgD) of the population) as the following formulas:(40)Vik+1=w·Vik+c1·rand·Pi-Xik+c2·rand·Pg-Xik,Xik+1=Xik+Vik+1,where i=1,2,…,M, M is the swarm size, w is the inertia weight, c1 and c2 are positive constants called acceleration coefficients for regulating the maximum step size for the particle i to fly towards Pi and Pg, respectively, superscript k represents the algorithm iteration count, and rand() is a random value in the range [0,1]. Xi and Vi should be clamped to the range of [Xmin,Xmax] and [Vmin,Vmax] to control excessive roaming of particles outside the search space, respectively.Because of its ergodicity and stochasticity, the piecewise linear chaotic map is used to improve the global search capability of PSO, which is defined as follows:(41)cxik+1=cxikq,cxik∈0,q,1-cxik1-q,cxik∈q,1,where cxik+1 denotes the chaotic mapping of ith chaotic variable in the kth iteration and q∈(0,0.5)∪(0.5,1) is a given parameter. Under this chaotic mapping, cx behaves chaotically in (0,1). Moreover, based on the fact that chaotic search is efficient in small range, we take chaotic mapping within a dynamically shrinking range around the global best position of the swarm by the following formula:(42)x^idk+1=pgd·1+η·cxik+1,where d=1,2,…,D, x^idk+1 denotes the chaotic mapping of the particle i in the dth dimension, pgd is the global best position in the dth dimension, and η is a scaling parameter with an empirical value 1.1. After chaos, the particle i should still be clamped to the range of [Xmin,Xmax] by the following formula:(43)Xik+1=max⁡min⁡X^ik+1,Xmax,Xmin.Integrated with chaotic map and SQP, the hybrid optimization algorithm CPSO-SQP is described as follows.Step 1.
Initialize swarm size, maximum function evaluations, other parameters, and the initial positionXi and velocity Vi of each particle in the swarm.Step 2.
Evaluate the objective function of each particle, then set the initial personal best solutionPi, select the global best particle Pg of CPSO to guide the swarm evolution direction, and set Po=Pg which is the global best particle of CPSO-SQP.Step 3.
Update the velocity and position of each particle according to (40) if the particle is not in stagnation, otherwise, according to (41)–(43).Step 4.
Evaluate the objective function of each particle in the new position.Step 5.
If the current position is better than personal best solution, updatePi; otherwise, set this particle in stagnation.Step 6.
IfPi is better than global best solution of CPSO, update Pg and set a flag for arousing SQP. Moreover, if Pg is better than Po, update Po.Step 7.
IfPg is just updated (the flag for SQP is true), set Pg as initial value and run SQP. If the result of SQP is better than Po, update Po.Step 8.
If maximum function evaluation is reached, outputPo as the global best solution; otherwise, go back to Step 3 for continuing iterations.In this algorithm, the powerful global search capability of CPSO to search the solution space and the solution obtained from CPSO is taken as an initial estimate of SQP algorithm to compensate for the weakness of SQP weak global search. With good local convergence performance and strong nonlinear convergence speed, SQP algorithm makes a careful search for the optimization problem to make up shortcomings of the CPSO weak local search. Thus, the possibility of exploring a global minimum in problems with more local optima is increased. Benefiting from the fast globally converging characteristics of CPSO and the effective local search ability of SQP, the proposed method can obtain the global optimal result of our problem through constant iteration.


## 5. Numerical Simulation

In this section, we present some numerical simulation results by the proposed method in this paper. Takea=b=γ=k=1, l=1, and r=2. Moreover, let (44)fx,t=e-tcos⁡x,gx,t=-π2sin⁡πtcos⁡πx,w0x=sin⁡πx,w1x=0,θ0x=sin⁡x.Thus, if we take U(t)=(u1(t),u2(t))⊤=0 for t>0, the uncontrolled equation (1) can be solved analytically and the exact state solution is (45)wx,t=sin⁡πxcos⁡πt,θ=e-tsin⁡x.Obviously, direct calculation yields that (46)E0=w0,wt0,θ0⊤H=π22+e-2t12-14sin⁡21/2t=0=2.28>rwhich implies that the constraint of energy function is not satisfied at the initial time.In this section, we take the admissible control setU as the form(47)U=U·=u1,u2T∈L∞0,∞×L∞0,∞∣-a≤uit≤aa.e.t∈0,∞,i=1,2,where positive constant a is the bound of the control constraint.For the spatial discretization of the thermoelastic system, we use the basis functionsϕi(x) for any 1≤i≤N as follows: (48)ϕix=Nx-xi-1,xi-1≤x<xi,-Nx-xi+1,xi≤x<xi+1,0,otherwise,where N is the number of subintervals in the spatial domain. Here, we choose N=20, which means that the finite element triangulation Th divides the space domain [0,l] into 20 subintervals.Through optimal time parameterization, the problem(TOCPh) has a fixed time horizon [0,1]. Using control parameterization, we subdivide the time interval [0,1] into p subintervals with p=20 and approximate the control functions by the piecewise-constant function in each subinterval [tk-1,tk),k=1,2,…,p. Then, the derived optimal parameter selection problem can be solved by the CPSO-SQP algorithm.Our numerical simulation study was carried out within the MATLAB programming environment running on a personal computer with the following configuration: Intel Core i5-3470 3.20 GHz CPU, 4 GB RAM, 64-bit Windows 7 Operating System.Figure1 shows the relationship of the optimal time T∗ and the control bound a. As the increase of a, the shortest time T∗, which is needed to drive the solution of the system to the energy constrained range, drops clearly. This means that the control ability of u∗ with large a is stronger than the that with small a. In Figure 2, the computed optimal control is presented for a=1,2,3,4, which shows that the Bang-Bang property holds for these cases. Moreover, we still make the comparison of states w and θ with control u=0 and the optimal states w∗ and θ∗ with control u∗ for the case a=4 in Figure 3.Figure 1
The optimal timeT∗ for different control bound a.Figure 2
The comparison of optimal controls fora=1,2,3,4.
(a)
a
=
1
(b)
a
=
2
(c)
a
=
3
(d)
a
=
4Figure 3
The comparison of statesw and θ with control u=0 and states w∗ and θ∗ with control u∗ for a=4.
(a)
Uncontrolled initial statew
(b)
Uncontrolled initial stateθ
(c)
Controlled statew
(d)
Controlled stateθ


## 6. Conclusions

In this paper, we have proposed a numerical approximation method for the time optimal control problem of the thermoelastic system. Following the schedules of Galerkin finite element discretization, time parameterization, and control parameterization, the original optimal control problem was reduced to an optimal parameter selection problem governed by a system of ordinary differential equations. In this problem, the optimal time and corresponding optimal controls are taken as decision variable to be optimized. Finally, a hybrid optimization algorithm CPSO-SQP was used to solve this problem. Simulation results demonstrate that the proposed numerical approximation method is effective at driving the thermoelastic system state to a target domain in the shortest time. Nevertheless, there is still room for improvement. We used fixed knot points for the control discretization in Section4. For more accurate results, it is possible to take the control switching points as decision variables by the so-calledtime-scaling transformation described in [38].



---

*Source: 6047670-2017-05-24.xml*