A Numerical Method of the Euler-Bernoulli Beam with Optimal Local Kelvin-Voigt Damping

Abstract

This paper deals with the numerical approximation problem of the optimal control problem governed by the Euler-Bernoulli beam equation with local Kelvin-Voigt damping, which is a nonlinear coefficient control problem with control constraints. The goal of this problem is to design a control input numerically, which is the damping and distributes locally on a subinterval of the region occupied by the beam, such that the total energy of the beam and the control on a given time period is minimal. We firstly use the finite element method (FEM) to obtain a finite-dimensional model based on the original PDE system. Then, using the control parameterization method, we approximate the finite-dimensional problem by a standard optimal parameter selection problem, which is a suboptimal problem and can be solved numerically by nonlinear mathematical programming algorithm. At last, some simulation studies will be presented by the proposed numerical approximation method in this paper, where the damping controls act on different locations of the Euler-Bernoulli beam.

1. Introduction

Let 𝑇,𝐿>0 be two positive constants. We denote by 𝑄 the product set (0,𝐿)×(0,𝑇). Consider a nonhomogeneous clamped elastic beam of length 𝐿, where one segment of the beam is made of a viscoelastic material with Kelvin-Voigt constitutive relation. By the Kirchhoff hypothesis, neglecting the rotatory inertia, the transversal vibration of the beam can be described by the following equation and boundary-initial conditions: 𝜌𝑤̈+𝑞𝑤^′′+𝐷𝑤̇^′′′′()𝑤()()=𝑓,𝑥,𝑡∈𝑄,0,𝑡=𝑤𝐿,𝑡=𝑤^′(0,𝑡)=𝑤^′()𝑤()𝐿,𝑡=0,𝑥,0=𝑤₀(𝑥),𝑤̇(𝑥,0)=𝑤₁(𝑥), (1) where 𝑤 represents transversal displacement of the beam, 𝑤₀,𝑤₁∈𝐿²(0,𝐿) are given initial data, and the notations 𝑤̇ and 𝑤^′ denote the derivatives with respect to the temporal variable and the spatial variable of 𝑤, respectively. Here 𝑓∈𝐿²(𝑄) is an external applied distributed force, 𝜌 is the linear mass density of the beam material, 𝑞 is the flexural rigidity, and 𝐷 is the Kelvin-Voigt damping coefficient. In this paper, we assume 𝜌,𝑞∈𝐿^∞(0,𝐿), 𝐷∈𝐿^∞(𝑄), and 𝜌,𝑞≥𝑐>0 for 𝑥∈(0,𝐿) for some constant 𝑐. we say that the Kelvin-Voigt damping is globally distributed if the damping coefficient 𝐷≥𝑐>0 on [0,𝐿]; we say it is locally distributed if 𝐷≥𝑐>0 only on some subinterval of [0,𝐿] and 𝐷=0 elsewhere.

Smart materials such as shape memory alloys and piezoceramics [1–3] have been applied for active vibration control of elastic structures. Accordingly, one can introduce control terms to the elastic systems such as the damping coefficients and Young’s moduli as well. Magnetorheological (MR) dampers [4–7] are one of the most promising new actuation mechanisms that use MR fluids to provide variable damping actuation for active control of structures. Because of their mechanical simplicity, high dynamic range, low power requirements, large force capacity, and robustness, these devices have been shown to mesh well with application demands and constraints to offer an attractive control method to structural vibration.

In this paper, we will study the Euler-Bernoulli beam equation with optimal local Kelvin-Voigt damping. To be more specific, let 𝜔∈[0,𝐿] be a subinterval and let 𝜒_𝜔 be the characteristic function of 𝜔, that is, 𝜒_𝜔(𝑥)=1,𝑥∈𝜔,0,otherwise; (2) We define 𝒰=𝑢(⋅)∈𝐿^∞(0,𝑇)∣𝑅₁≤𝑢(𝑡)≤𝑅₂fora.e.𝑡∈[0,𝑇] (3) and assume the damping coefficient 𝐷 has the following form: 𝐷(𝑥,𝑡)=𝜒_𝜔(𝑥)𝑢(𝑡),(𝑥,𝑡)∈𝑄, (4) where 𝑅₁,𝑅₂∈ℝ are two fixed constants and 𝑢∈𝒰 is a control function.

Let 𝑉=𝐻²₀(0,𝐿) with the norm ‖𝑤‖_𝑉=^𝐿₀𝑞(𝑥)𝑤^′′(𝑥)²𝑑𝑥^1/2,∀𝑤∈𝑉, (5) and 𝐻=𝐿²_𝜌(0,𝐿) with the norm ‖𝑣‖_𝐻=^𝐿₀𝜌(𝑥)|𝑣(𝑥)|²𝑑𝑥^1/2,∀𝑣∈𝐻. (6) Define ℋ=𝑉×𝐻 with the norm 𝑤𝑣_ℋ=‖𝑤‖²_𝑉+‖𝑣‖²_𝐻^1/2𝑤𝑣,∀∈ℋ. (7) Then, ℋ is Hilbert space and the energy of the beam at time 𝑡 is 𝑤(𝑡)𝑤̇(𝑡)²_ℋ=^𝐿₀𝑞(𝑥)𝑤^′′(𝑥,𝑡)²+𝑑𝑥_𝐿0𝜌(𝑥)|𝑤̇(𝑥,𝑡)|²𝑑𝑥, (8) where 𝑤 is the solution of (1). The optimal control problem that we will study is formulated as follows: (OCP)min_𝑢∈𝒰𝐽(𝑢)=12^𝑇₀𝑤(𝑡)𝑤̇(𝑡)²_ℋ+1𝑑𝑡2^𝑇₀|𝑢(𝑡)|²𝑑𝑡 (9) subject to the controlled equation (1), where 𝑤 is the solution of (1). Throughout the paper, we will omit the notations 𝑡 or 𝑥 in the functions of 𝑡 or 𝑥 in the case that there is no risk to make any confusion.

Due to the importance from both perspectives of mathematics and applied science, the control problem of various beam equations has been considered by many researchers [8–10]. The study of the Euler-Bernoulli beam is one of the most active research topics in control theory. In [11], the authors consider the vibration of the Euler-Bernoulli beam with Kelvin-Voigt damping distributed locally on any subinterval of the region occupied by the beam. By making use of the frequency domain method and the multiplier technique, they prove that the semigroup associated with the equation for the transversal motion of the beam is exponentially stable. In [12], the author studies the basis property and the stability of a distributed system described by a nonuniform Euler-Bernoulli beam equation under linear boundary feedback control. The Riesz basis property is presented and the exponential stability is concluded. In [13], stabilization of Euler-Bernoulli beam by means of the pointwise feedback force is considered. Both uniform and nonuniform energy decay may occur, which depend on the boundary conditions. There are some other related papers about the studies of Euler-Bernoulli beam equations [14–17].

In this paper, we will study the numerical approximation of the optimal control problem (OCP), which is a nonlinear bilinear control problem with control constraint. Bilinear control problems are already studied by many researchers [18–20]. In our paper, we want to design a damping control numerically, which acts on local interval of the beam, such that the total energy of the beam and the control on a given time period is minimal. It appears that little work has been done on numerical methods for this problem. By the standard finite element method (FEM), problem (OCP) was firstly approximated by an optimal control problem (OCP^ℎ) governed by a system of ordinary differential equations. Then, using the control parameterization method [21], we will approximate the finite-dimensional problem by another suboptimal problem (OCP^ℎ_𝑑), which is a standard optimal parameter selection problem and can be solved numerically by nonlinear mathematical programming algorithm. At last, some simulation studies will be presented by the numerical method proposed in this paper.

2. The Semidiscrete Approximation by FEM

In this section, we will approximate the original optimal control problem (OCP) with FEM method. Noting that (1) involves the spatial derivative of four orders, the conforming FEM space should belong to 𝐻²(0,1). Consider the interval domain [0,𝐿]. The triangulation 𝒯^ℎ of [0,𝐿] divides [0,𝐿] into a finite number of subintervals 𝐼_𝑗=[𝑥_𝑗−1,𝑥_𝑗], 𝑗=1,2,…,𝑁+1, using the grid points: 0=𝑥₀<𝑥₁<⋯<𝑥_𝑁<𝑥_𝑁+1=𝐿, (10) where we will call 𝐼_𝑗=[𝑥_𝑗−1,𝑥_𝑗] the 𝑗th element and ℎ_𝑗=𝑥_𝑗−𝑥_𝑗−1 the size of this element. The discretization parameter ℎ is the maximum size of all ℎ_𝑗, 𝑗=1,2,…,𝑁+1. Associated with every triangulation 𝒯^ℎ, we define a finite-dimensional space as follows: 𝑉^ℎ=𝑣^ℎ∈𝐶¹[]∣𝑣0,𝐿_ℎ∈𝑃₃𝐼_𝑗𝑣for𝑗=1,2,…,𝑁+1,_ℎ(0)=𝑣^ℎ(𝐿)=𝑣^ℎ′(0)=𝑣^ℎ′(𝐿)=0, (11) where 𝑃₃(𝐼_𝑗) is the space of all polynomials of degree less than or equal to 3 over the subinterval 𝐼_𝑗. Obviously, we have 𝑉^ℎ∈𝐻²₀(0,𝐿). Thus, we can write 𝑉^ℎ=span𝜙_𝑖,∣𝑖=1,2,…,2𝑁, (12) where, for 𝑖=1,2,…,𝑁, 𝜙_𝑖𝑥_𝑘=𝛿_𝑖𝑘,𝜙^′_𝑖𝑥_𝑘=0,∀𝑘=0,1,…,𝑁+1, (13) and, for 𝑖=𝑁+1,𝑁+2,…,2𝑁, 𝜙_𝑖𝑥_𝑘=0,𝜙^′_𝑖𝑥_𝑘=𝛿_{𝑖−𝑁,𝑘},∀𝑘=0,1,…,𝑁+1. (14)

Define a bilinear form (⋅,⋅)_𝜌 over 𝐿²_𝜌(0,𝐿)×𝐿²_𝜌(0,𝐿) by setting 𝑓,𝑔_𝜌=^𝐿₀𝜌(𝑥)𝑓(𝑥)𝑔(𝑥)𝑑𝑥,∀𝑓,𝑔∈𝐿²_𝜌(0,𝐿). (15) Define another two bilinear forms 𝑎_𝑞(⋅,⋅) and 𝑎_𝜔(⋅,⋅) over 𝐻²₀(Ω)×𝐻²₀(Ω) by setting 𝑎_𝑞𝑓,𝑔=^𝐿₀𝑞(𝑥)𝑓^′′(𝑥)𝑔^′′(𝑥)𝑑𝑥,∀𝑓,𝑔∈𝐻²₀(Ω),𝑎_𝜔𝑓,𝑔=_𝜔𝑓^′′(𝑥)𝑔^′′(𝑥)𝑑𝑥,∀𝑓,𝑔∈𝐻²₀(Ω), (16)(17) respectively. Obviously, the two bilinear forms (⋅,⋅)_𝜌 and 𝑎_𝑞(⋅,⋅) are the inner products of 𝐿²_𝜌(0,𝐿) and 𝐻²₀(Ω), respectively. Then, the finite element approximation of (1) consists in finding 𝑤^ℎ(𝑡)=𝑤^ℎ(⋅,𝑡), which belongs to 𝑉^ℎ for 𝑡∈[0,𝑇], and satisfies 𝑤̈^ℎ(𝑡),𝑣^ℎ_𝜌+𝑎_𝑞𝑤^ℎ(𝑡),𝑣^ℎ+𝑢(𝑡)𝑎_𝜔𝑤̇^ℎ,𝑣^ℎ=𝑓(𝑡),𝑣^ℎ∀𝑣_ℎ∈𝑉^ℎ𝑤,0<𝑡≤𝑇,_ℎ(0)=𝑤^ℎ₀,𝑤̇^ℎ(0)=𝑤^ℎ₁, (18) where (⋅,⋅) denotes the standard inner product of 𝐿²(0,𝐿) and 𝑤^ℎ₀,𝑤^ℎ₁∈𝑉^ℎ are the proper approximations of 𝑤₀,𝑤₁ on 𝑉^ℎ. In the following, we write 𝑤^ℎ(𝑥,𝑡)=^2𝑁_𝑗=1𝑋_𝑗(𝑡)𝜙_𝑗(𝑥),𝑤^ℎ(𝑥,0)=𝑤^ℎ₀(𝑥)=^2𝑁_𝑗=1𝑋_0𝑗𝜙_𝑗(𝑥),𝑤̇^ℎ(𝑥,0)=𝑤^ℎ₁(𝑥)=^2𝑁_𝑗=1𝑌_0𝑗𝜙_𝑗(𝑥). (19)(20) Substituting (19) into (18) and taking 𝑣_ℎ=𝜙_𝑖 yield that ^2𝑁_𝑗=1𝜙_𝑖,𝜙_𝑗𝜌𝑋̈_𝑗(𝑡)+^2𝑁_𝑗=1𝑎_𝑞𝜙_𝑖,𝜙_𝑗𝑋_𝑗(𝑡)+𝑢(𝑡)^2𝑁_𝑗=1𝑎_𝜔𝜙_𝑖,𝜙_𝑗𝑋̇_𝑗(𝑡)=𝑓(𝑡),𝜙_𝑖. (21) Moreover, in this paper, we take 𝑤^ℎ₀,𝑤^ℎ₁ as the 𝐿²_𝜌-projection approximations of 𝑤₀,𝑤₁ on 𝑉^ℎ; that is, 𝑤^ℎ₀,𝑣^ℎ_𝜌=𝑤₀,𝑣^ℎ_𝜌,𝑤^ℎ₁,𝑣^ℎ_𝜌=𝑤₁,𝑣^ℎ_𝜌,∀𝑣^ℎ∈𝑉^ℎ. (22) Then, substituting (20) into (22) and taking 𝑣_ℎ=𝜙_𝑖 yield that ^2𝑁_𝑗=1𝜙_𝑖,𝜙_𝑗𝜌𝑋_0𝑗=𝑤₀,𝜙_𝑖𝜌,^2𝑁_𝑗=1𝜙_𝑖,𝜙_𝑗𝜌𝑌_0𝑗=𝑤₁,𝜙_𝑖𝜌. (23) Define 𝑋(𝑡)=𝑋_𝑖(𝑡)_2𝑁×1,𝑀=𝜙_𝑖,𝜙_{𝑗𝜌2𝑁×2𝑁},𝐾=𝑎_𝑞𝜙_𝑖,𝜙_{𝑗2𝑁×2𝑁},𝑅=𝑎_𝜔𝜙_𝑖,𝜙_{𝑗2𝑁×2𝑁},𝑋₀=𝑋_{0𝑖2𝑁×1},𝑌₀=𝑌_{0𝑖2𝑁×1},𝐹(𝑡)=𝑓(𝑡),𝜙_𝑖2𝑁×1,𝜂=𝑤₀,𝜙_{𝑖𝜌2𝑁×1},𝜉=𝑤₁,𝜙_{𝑖𝜌2𝑁×1}. (24) Thus, by (21) and (23), we can obtain the following system of controlled ordinary differential equations: 𝑀𝑋̈(𝑡)(𝑡)(𝑡)𝑅𝑋̇(𝑡)(𝑡),𝑋(0)=−𝐾𝑋−𝑢+𝐹=𝑋₀,𝑋̇(0)=𝑌₀, (25) where 𝑋₀=𝑀⁻¹𝜂, 𝑌₀=𝑀⁻¹𝜉. Let ̇𝑋𝑌𝑋=𝑌,𝑍=. (26) We define 𝐺(𝑡,𝑍(𝑡),𝑢(𝑡))=0𝐼−𝑀⁻¹𝐾−𝑢(𝑡)𝑀⁻¹𝑅𝑍(𝑡)+0𝑀⁻¹𝐹(𝑡). (27) Then (25) can be rewritten as 𝑍̇(𝑡)((𝑡)(𝑡,𝑍(0)=𝐺𝑡,𝑍,𝑢))=𝑍₀, (28) where 𝑍₀=^𝑋0_𝑌0.

By (19), a direct computation yields 𝑤^ℎ(𝑡)²_𝑉=𝑎_𝑞𝑤^ℎ(𝑡),𝑤^ℎ(𝑡)=^2𝑁_𝑖=1^2𝑁_𝑗=1𝑎_𝑞𝜙_𝑖,𝜙_𝑗𝑋_𝑖(𝑡)𝑋_𝑗(𝑡)=(𝑋,𝐾𝑋)_ℝ^2𝑁, (29) where (⋅,⋅)_ℝ^2𝑁 denotes the inner product of ℝ^2𝑁 and 𝑉=𝐻²₀(0,𝐿). Similarly, we have 𝑤̇^ℎ(𝑡)²_𝐻=𝑤̇^ℎ(𝑡),𝑤̇^ℎ(𝑡)_𝜌=^2𝑁_𝑖=1^2𝑁_𝑗=1𝜙_𝑖,𝜙_𝑗𝜌𝑋̇_𝑖(𝑡)𝑋̇_𝑗(𝑡)=(𝑌,𝑀𝑌)_ℝ^2𝑁, (30) where 𝐻=𝐿²_𝜌(0,𝐿). Define Π=.𝐾00𝑀 (31) Then, it follows from (29), (30), and (31) that 𝐽(𝑢)≈𝐽^ℎ(𝑢)=12^𝑇₀𝑤^ℎ(𝑡)𝑤̇^ℎ(𝑡)²_ℋ1𝑑𝑡+2^𝑇₀|𝑢(𝑡)|²=1𝑑𝑡2^𝑇₀(𝑍(𝑡),Π𝑍(𝑡))_ℝ^4𝑁1𝑑𝑡+2^𝑇₀|𝑢(𝑡)|²𝑑𝑡. (32) Thus, by combining (28) and (32), the semidiscrete approximation of problem (OCP) is formulated as follows: OCP^ℎmin_𝑢∈𝒰𝐽^ℎ(𝑢) (33) subject to (28).

3. Piecewise-Constant Control Approximation

In general, problem (OCP^ℎ) cannot be solved analytically. Using the control parameterization method, which has been successfully applied to provide numerical solutions for a wide variety of practical optimal control problems [21–24], we will approximate problem (OCP^ℎ) by a standard optimal parameter selection problem. This method involves approximating the control function by a piecewise-constant function with possible discontinuities at a set of preassigned switching points, which produces an approximation problem such that the solution of this approximation is a suboptimal solution to problem (OCP^ℎ).

Let 𝜏_𝑘, 𝑘=0,1,…,𝑑, be prefixed time knot points satisfying 0=𝜏₀<𝜏₁<𝜏₂<⋯<𝜏_𝑑=𝑇. (34) With piecewise-constant basis functions, the control input 𝑢 for the problem (OCP^ℎ) is approximated over the 𝑘th control subinterval [𝜏_𝑘−1,𝜏_𝑘) as follows: 𝑢(𝑡)≈𝑢_𝑑(𝑡)=𝜇_𝑘,𝑡∈𝜏_𝑘−1,𝜏_𝑘,𝑘=1,2,…,𝑑, (35) where 𝜇_𝑘 is the value of the control on the 𝑘th subinterval [𝜏_𝑘−1,𝜏_𝑘). Define 𝒰_𝑑=𝜇=𝜇₁,𝜇₂,…,𝜇_𝑑^⊤∈ℝ^𝑑∣𝑅₁≤𝜇_𝑘≤𝑅₂,𝑘=1,2,…,𝑑. (36) Then the approximate piecewise-constant control can be written as follows: 𝑢_𝑑𝑡∣𝜇=^𝑑_𝑘=1𝜇_𝑘𝜒_{[𝜏𝑘−1,𝜏𝑘)}(𝑡),𝑡∈[0,𝑇],𝜇∈𝒰_𝑑, (37) where 𝜒_{[𝜏𝑘−1,𝜏𝑘)} is the characteristic function of the interval [𝜏_𝑘−1,𝜏_𝑘), 𝑘=1,2,…,𝑑. Substituting (37) into the dynamic system (28) yields that 𝑍̇(𝑡)=𝐺𝑡,𝑍(𝑡),𝑢_𝑑=𝑡∣𝜇_𝑑𝑘=1𝐺𝑡,𝑍(𝑡),𝜇_𝑘𝜒_{[𝜏𝑘−1,𝜏𝑘)}(𝑡)[],𝑍(0),𝑡∈0,𝑇=𝑍₀. (38) Let 𝑍_𝑑(⋅∣𝜇) denote the solution of system (38) corresponding to 𝜇∈𝒰_𝑑. Thus, from the problem (OCP^ℎ), we can obtain another parameter optimization problem, which is stated as follows: OCP^ℎ_𝑑min_{𝑢∈𝒰𝑑}𝐽^ℎ_𝑑𝜇 (39) subject to (38), where 𝐽^ℎ_𝑑𝜇=12^𝑇₀𝑍_𝑑𝑡∣𝜇,Π𝑍_𝑑𝑡∣𝜇_ℝ4𝑁+1𝑑𝑡2^𝑇₀𝑢_𝑑𝑡∣𝜇²𝑑𝑡. (40)

After the parameterization of control, problem (OCP^ℎ_𝑑) involves a finite number of decision variables. Thus, it should be much easier to solve than problem (OCP^ℎ), which involves determining the value of a function at an infinite number of time points.

4. Variational Method for Solving Problem (OCP^ℎ_𝑑)

Problem (OCP^ℎ_𝑑) is an optimal parameter selection problem in the canonical form [24], which can be solved as nonlinear optimization problems using the SQP method. Standard SQP algorithm for nonlinear optimization exploits the gradient of the cost functional to generate search directions that lead to profitable areas of the search space [25, 26]. For the approximate problem (OCP^ℎ_𝑑), the cost functional is implicit function of the decision vector 𝜇. Using the variational method [27, 28], we can compute this gradient and solve problem (OCP^ℎ_𝑑).

For each 𝑚=1,2,…,𝑑, it follows from (38) that 𝑍_𝑑𝑡∣𝜇=𝑍_𝑑𝜏_𝑚−1+∣𝜇_{𝑡𝜏𝑚−1}𝐺𝑡,𝑍_𝑑𝑠∣𝜇,𝜇_𝑚𝜏𝑑𝑠,𝑡∈_𝑚−1,𝜏_𝑚. (41) Then for 𝑡∈[𝜏_𝑚−1,𝜏_𝑚), differentiating (41) with respect to 𝜇_𝑘 yields that 𝜕𝑍_𝑑𝑡∣𝜇𝜕𝜇_𝑘=𝜕𝑍_𝑑𝜏_𝑚−1∣𝜇𝜕𝜇_𝑘+^𝑡_𝜏𝑚−1𝜕𝐺𝑠,𝑍_𝑑𝑠∣𝜇,𝜇_𝑚𝜕𝑍_𝑑𝜕𝑍_𝑑𝑠∣𝜇𝜕𝜇_𝑘𝑑𝑠,for𝑘<𝑚,𝜕𝑍_𝑑𝑡∣𝜇𝜕𝜇_𝑘=𝜕𝑍_𝑑𝜏_𝑚−1∣𝜇𝜕𝜇_𝑘+^𝑡_𝜏𝑚−1𝜕𝐺𝑠,𝑍_𝑑𝑠∣𝜇,𝜇_𝑚𝜕𝑍_𝑑𝜕𝑍_𝑑𝑠∣𝜇𝜕𝜇_𝑘+𝑑𝑠_{𝑡𝜏𝑚−1}𝜕𝐺𝑠,𝑍_𝑑𝑠∣𝜇,𝜇_𝑚𝜕𝑢𝑑𝑠,for𝑘=𝑚,𝜕𝑍_𝑑𝑡∣𝜇𝜕𝜇_𝑘=0,for𝑘>𝑚. (42) Here we have 𝜕𝐺𝑡,𝑍_𝑑𝑡∣𝜇,𝜇_𝑚𝜕𝑍_𝑑=0𝐼−𝑀⁻¹𝐾−𝜇_𝑚𝑀⁻¹𝑅,𝜕𝐺𝑡,𝑍_𝑑𝑡∣𝜇,𝜇_𝑚=𝜕𝑢000−𝑀⁻¹𝑅𝑍_𝑑𝑡∣𝜇, (43) for 𝑡∈[𝜏_𝑚−1,𝜏_𝑚). Define 𝛿_𝑘𝑚=𝛿1,if𝑘=𝑚,0,otherwise,_𝑘𝑚=1,if𝑘≤𝑚,0,otherwise. (44) Then, by (42) we obtain 𝜕𝑍_𝑑𝑡∣𝜇𝜕𝜇_𝑘=𝛿_𝑘𝑚𝜕𝑍_𝑑𝜏_𝑚−1∣𝜇𝜕𝜇_𝑘+𝛿_𝑘𝑚^𝑡_𝜏𝑚−1𝜕𝐺𝑠,𝑍_𝑑𝑠∣𝜇,𝜇_𝑚𝜕𝑍_𝑑𝜕𝑍_𝑑𝑠∣𝜇𝜕𝜇_𝑘𝑑𝑠+𝛿_𝑘𝑚^𝑡_𝜏𝑚−1𝜕𝐺𝑠,𝑍_𝑑𝑠∣𝜇,𝜇_𝑚𝜏𝜕𝑢𝑑𝑠,𝑡∈_𝑚−1,𝜏_𝑚. (45) Moreover, it is easy to see that 𝜕𝑍_𝑑0∣𝜇𝜕𝜇_𝑘=0. (46) For 𝑡∈[𝜏_𝑚−1,𝜏_𝑚), differentiating (45) with respect to 𝑡 yields that 𝑑𝑑𝑡𝜕𝑍_𝑑𝑡∣𝜇𝜕𝜇_𝑘=𝛿_𝑘𝑚𝜕𝐺𝑡,𝑍_𝑑𝑡∣𝜇,𝜇_𝑚𝜕𝑍_𝑑𝜕𝑍_𝑑𝑡∣𝜇𝜕𝜇_𝑘+𝛿_𝑘𝑚𝜕𝐺𝑡,𝑍_𝑑𝑡∣𝜇,𝜇_𝑚.𝜕𝑢 (47)

Now, we define Γ_𝑘𝑡∣𝜇=𝜕𝑍_𝑑𝑡∣𝜇𝜕𝜇_𝑘,𝑘=1,2,…,𝑑. (48) Then, it follows from (46) and (47) that 𝑑Γ_𝑘𝑡∣𝜇=𝛿𝑑𝑡^𝑘𝑚𝜕𝐺𝑡,𝑍_𝑑𝑡∣𝜇,𝜇_𝑚𝜕𝑍_𝑑Γ_𝑘𝑡∣𝜇+𝛿_𝑘𝑚𝜕𝐺𝑡,𝑍_𝑑𝑡∣𝜇,𝜇_𝑚,𝜏𝜕𝑢𝑡∈_𝑚−1,𝜏_𝑚,𝑚=1,2,…,𝑑, (49) with the initial condition Γ_𝑘0∣𝜇=0. (50) As a result, by using the chain rule, we can derive the gradient of 𝐽^ℎ_𝑑(𝜇) with respect to 𝜇_𝑘, 𝑘=1,2,…,𝑑, as follows: 𝜕𝐽^ℎ_𝑑𝜇𝜕𝜇_𝑘=^𝑇₀𝑍_𝑑𝑡∣𝜇,ΠΓ_𝑘𝑡∣𝜇_ℝ4𝑁𝑑𝑡+𝜇_𝑘𝜏_𝑘−𝜏_𝑘−1. (51) By incorporating these formulae into the SQP algorithm, we can solve the problem (OCP^ℎ_𝑑) numerically.

5. Numerical Simulations

In this section, we present some numerical simulation results by the approximation method presented in this paper. Let 𝐿=𝑇=1, 𝜌(𝑥)=𝑞(𝑥)=1, 𝑅₁=0, and 𝑅₂=1. Moreover, we take 𝑤₀(𝑥)=10𝑥²(1−𝑥)²,𝑤₁(𝑥)𝑓()==0,𝑥,𝑡240−10𝜋²𝑥²(1−𝑥)²cos(𝜋𝑡). (52) If the damping coefficient 𝐷(𝑥,𝑡)=0,(𝑥,𝑡)∈𝑄, (53) which implies that the system (1) is without any damping, the exact solution (see Figure 1) of (1) is 𝑤(𝑥,𝑡)=10cos(𝜋𝑡)𝑥²(1−𝑥)²,(𝑥,𝑡)∈𝑄. (54) In this case, we have 12^𝑇₀𝑤(𝑡)𝑤̇(𝑡)²_ℋ1𝑑𝑡+2¹₀|𝑢(𝑡)|²𝑑𝑡=20.39. (55)

Figure 1 

The state without damping control; that is, 𝑢=0.

In the following, we discuss the optimal control problem (OCP) with local Kelvin-Voigt damping acting on two different locations. We write 𝜔₁=[0.6,1] and 𝜔₂=[0.2,0.6]. By the approximation method presented in this paper, we can get the numerical results for the optimal control problem (OCP) with 𝜔=𝜔₁ and 𝜔=𝜔₂, which will be called problem (OCP)₁ and problem (OCP)₂, respectively. The numerical optimal controls for the two problems are presented by Figure 2. Moreover, the difference of 𝑤₁ and 𝑤 is shown in Figure 3, where 𝑤 is defined by (54) and 𝑤₁ is the numerical state function of (OCP)₁. The difference of 𝑤₂ and 𝑤 has the same property, where 𝑤₂ is the numerical state function of (OCP)₂. For problem (OCP)₁, the optimal value is 14.22 and for problem (OCP)₂, the optimal value is 15.65, which means that the control effect of the location 𝜔₁ is better than the location 𝜔₂.

Figure 2 

The optimal controls for (OCP)1 and (OCP)2.

Figure 3 

The difference between 𝑤1 and 𝑤.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China (61374096 and 61104048) and the Natural Science Foundation of Zhejiang (Y6110751).