Finite Element Approximation of Optimal Control Problem with Weighted Extended B-Splines

: In this research article, an optimal control problem (OCP) with boundary observations is approximated using ﬁnite element method (FEM) with weighted extended B-splines (WEB-splines) as basis functions. This type of OCP has a distinct aspect that the boundary observations are outward normal derivatives of state variables, which decrease the regularity of solution. A meshless FEM is proposed using WEB-splines, deﬁned on the usual grid over the domain, R 2 . The weighted extended B-spline method (WEB method) absorbs the regularity problem as the degree of the B-splines is increased. Convergence analysis is also performed by some numerical examples.


Introduction
The optimal control problems (OCPs) originate from the civil engineering of earth supplies [1]. OCPs governed by the partial differential equations (PDEs) feature in the practical applications of electrical, mechanical, control and aerospace engineering from many decades. Some large space structures, sensitive equipment and their operating structures also require optimal control technology. Although, civil engineering assembles the applications of OCPs in an exclusive style, which offers new challenges [2].
Since OCPs have widespread applications, so it is necessary to solve them by numerically stable methods whenever it is problematic to solve analytically. The numeric simulations of above stated applications have been performed using well known numeric techniques [3]. No doubt, practically, in order to avoid the unappropriate disturbances, closed loop optimal controls have a leading position over the open loop controls. Under some slight conditions on the coefficients and admissible control points, Sun et al. [4] gave the solution of closed loop OCP via Riccati equation.
It can easily be estimated from the increasing publishing trend of research articles over the last few decades that FEM has proved itself a well programmed and commonly used method in an extensive range of engineering applications in general, for example in structural analysis, heat transfer and fluid mechanics etc. Different mathematicians have presented FEM for PDEs as well as for OCPs [5][6][7][8][9][10][11][12].
Recently mathematicians used non uniform rational B-splines (NURBS) as the basis functions, instead of the customary Lagrange polynomials. Buckling analysis of thin symmetrically laminated composite plates using NURBS is done by [13]. Valizadeh et al. [14] studied the behaviour of functionally graded plates with NURBS as finite element basis functions. Liu et al. [15] analyzed the homogeneous and functionally graded microplates using an isogeometric analysis and non-classical Kirchhoff plate theory. The adaptive isogeometric analysis for two dimensional elasticity problems using locally refined B-splines is done by [16].
For optimal controls, Hinze et al. [17] studied the optimality conditions of OCPs in 2009. Becker et al. [18] analyzed constrained OCPs governed by the convection diffusion equations. FEM practiced the optimality conditions for approximation of problems governed by Stokes equations [19,20]. A variational discretization was done by [21]. Superconvergence properties of OCPs were firstly presented by [22]. As convergence is concerned, superconvergence for OCPs with semi-linear elliptic equations is done by [23]. The effort of [24] is also remarkable as it considered the elliptic OCPs with pointwise observations. Hou estimated errors of mixed FEM for elliptic OCPs [25]. Hou also discussed the error estimates of expanded mixed FEM for OCPs obtained by hyperbolic integro differential equation [26]. Chen et al. [27] applied mixed multiscale FEM for the convex OCPs with oscillating coefficients. Error estimates for mixed FEM for OCPs of low regularity are discussed in [28].
This manuscript presents FEM with WEB-splines to solve elliptic OCPs with boundary observations. The major characteristic of such optimal control problems is that the observations are the outward normal derivatives of the state variable on the boundary which decreases the regularity of solutions. There are certain limitations with standard FEM to deal with such type of inhomogeneous boundary conditions. In standard FEM, one is forced to approximate normal derivatives at boundary in the OCP by means of other numerical techniques, for example [29]. However mixed FEM solves such type of inhomogeneous boundary on a higher cost as compared to the standard FEM [29].
The WEB method resolves the regularity problem, relatively at a lower cost and boundary conditions are well approximated without being discontinuous. WEB method is a meshless technique defined on a usual grid over domain. Convergence rate of the WEB method depends on the degree of B-splines. A special extrapolation technique is employed to prevent the loss of stability in the presence of small cut-cell by the boundary. As a result, a new set of extended B-splines is obtained and is further stabilized by using the stabilizing factor. This stabilizing factor depends on the weight functions. Finally, the WEB-splines are obtained. For the error analysis of WEB method, a straight approach of finite element discretization could be tedious. The approximation order via quasi interpolant would be an easy approach to get the goal. Due to the standard projector property of canonical quasi interpolant, the task becomes easy because the approximation u h depends on Q q u within the support of B i , see (18). Approximation order of WEB method with elliptic boundary conditions gives the convergence of O(h d ), where d is the degree of B-spline [30]. But here we have an exceptional type of objective functional consisting of outward normal derivative of the state variable on boundary. Another main feature is that the adjoint state equation appears as an elliptic equation with inhomogeneous Dirichlet boundary condition whose value is the partial derivative of the state. Thus for an OCP of this type, WEB method can not approximate upto O(h d ). To the best of authors knowledge, such type of approximation to OCP with such special boundary observation has not been addressed earlier.
Let a square symmetric coefficient matrix A = a ij ∈ C (2,1) , defined onΩ be positive definite by nature. Then the OCP is to find the function y for which min u∈U ad Equation (2) is called the state equation, where domain Ω ⊂ R 2 is convex and has a smooth boundary ∂Ω. The fixed positive number β is a parameter of penalty. The function σ d ∈ L 2 (∂Ω) and ∂ n y is defined as A∇y · n where n is taken as outward drawn normal to boundary, f ∈ L 2 (Ω) and Ω ad is the admissible sub-domain of Ω.
is a continuous linear operator that maps the elements from L 2 (Ω ad ) to L 2 (Ω), which may be considered as a characteristic function χ Ω ad defined on Ω ad .
The set of admissible points U ad , is composed of admissible control points, u ∈ L 2 (Ω ad ) such that µ 1 ≤ u ≤ µ 2 , almost everywhere in Ω ad and µ 1 , µ 2 are real numbers.
Rest of the paper is organized in the following way. Section 2 contains some preliminary definitions. Also the OCP, optimality conditions and approximations are explained in this section. Section 3 represents the main hurdles in solving OCP by standard FEM and a detailed construction of WEB-splines over a bounded domain is presented in this section. Moreover, this section contains the pseudo-codes for application of WEB method over an OCP. Results regarding the approximation are given in Section 4. To support theoretical results, the numerical examples using proposed algorithm are given in Section 5. At last, some finishing conclusions are written in Section 6.

Preliminaries
Let Ω ⊂ R n , and we consider the continuous function z : Ω → R, making the functional space C(Ω). Let the multiindex α ∈ N n such that |α| = α 1 + ... + α n , and we define the space

For open and bounded domain
Classically, it is assumed that the functions in boundary value problems must have continuous derivatives. For PDEs, it is required to use Sobolev spaces W m,p (Ω) instead of the functional spaces C m (Ω). The typical representation of Sobolev spaces is clears the concept of weak derivatives where C ∞ c (Ω) is a set of continuous functions on a compact spaceΩ.
In Sobolev space, norm and semi norm are defined as z

this space
includes all of the functions in W m,p (Ω) that vanish on the boundary of Ω. Significantly, if we fix p = 2, Sobolev space W m,2 (Ω) can be represented as H m (Ω) and at the boundary W m,2 0 (Ω) can be represented as H m 0 (Ω). We denote the L 2 inner product defined on domain, boundary and admissible domain by (·, ·), ·, · and (·, ·) ad respectively.
The solution y of the boundary value problem lies in H 2 (Ω) ∩ H 1 0 (Ω) named as Y with normal derivative ∂ n y ∈ L 2 (∂Ω). The objective functional (1) is well defined.
The weak form of the control problem (1) is: find the solution y ∈ Y and u ∈ U ad satisfying the minimization problem min u∈U ad Here, a(·, ·) is the bilinear functional Y × H 1 0 −→ R and is symmetric as well as positive definite by nature.

Approximation of Model
The residual form of the problem (3) and (4) is represented as: Minimize I(y, u) subjects to the residual with Dirichlet boundary conditions y = 0 on ∂Ω. Consider the functions (y, u, σ) ∈ Y × U ad × L 2 (Ω), which are continuously differentiable while Y, U ad and L 2 (Ω) are Banach spaces. State equation r(y, u) = 0 carries a unique solution for all The following theorem shows the necessary and sufficient conditions for OCP. (3) and (4) of OCP exhibits the necessary and sufficient optimal conditions in Y × U ad × L 2 (Ω) stated as:

Theorem 1. The weak form Equations
with * is the adjoint of . Equivalently Equation (8) can be where P U ad is the orthogonal projection on to U ad .
The equivalent problems of the variational form of Equations (6)- (8) in the sense of distribution are the set of partial differential equations together with inequality constraint: Solving the set of Equations (13)-(15), we will have the solution of the control problem. It is easy to see that problem resembles to Dirichlet boundary control problem [29]. The state equation of the Dirichlet boundary control problem can be easily handled by using standard FEM as well as by the WEB method. The difficulty arises with the fact that boundary condition of the adjoint state equation is comprised of normal derivatives of the state equation. The WEB method addresses this problem effectively by using cubic WEB-splines.
The partial derivative of the state variable ∂ n y belongs to Bessel potential space lying on the boundary H 1/2 (∂Ω). Moreover, the norm of ∂ n y on H 1/2 (∂Ω) is bounded [29], i.e., The above relation shows the behavior of regularity of elliptic OCP with Dirichlet boundary conditions. The next lemma proves the regularity result of OCP and gives brief properties of functions involved in OCP.

Lemma 1.
Let the smooth bounded domain Ω, assuming f ∈ L 2 (Ω) and σ d ∈ H 1/2 (∂Ω) and let u ∈ L 2 (Ω ad ) then the solution of the state equation y and adjoint state equation σ has the properties: Proof. The properties of the Bessel potential space can easily prove the first property. For the second property, σ d ∈ H 1/2 (∂Ω) = θ o (H 1 (Ω)), where θ o is the trace operator which coincides with any usual restriction operator for continuous function. So, we may write σ d ∈ * (H 1 (Ω)). By using Equation (12), we get σ ∈ H 1 (Ω). Third property can be directly obtained from Equation (8).

WEB Based Finite Element Method
In standard FEM, a solution y is approximated by y h ∈ Y h by using variational formulation, where Y h is a finite dimensional linear space. The function y h , which is obtained from the system of algebraic linear equations, is continuous and piecewise linear on Ω h , where Ω h is some polygonal subdomain of Ω and vanishes on ∂Ω h . The triangular discretization is a very common method for the discretization of a domain. The triangulation is performed in such a way that for triangles λ i and λ j , i = j,λ i ∩λ j = Φ or none of the vertex of λ i lying in theλ j and λ i ∪ λ j , ∀ i, j establishes Ω h . Let h be the largest length of sides of triangles, built in the construction of Ω h . Triangulation becomes finer as h → 0 and the boundary ∂Ω h becomes closer to ∂Ω. To approximate OCP, state variable y is to be further approximated by y h , as its normal derivative in boundary condition of adjoint state equation does not allow the implementation of Dirichlet boundary condition directly. Yan et al. [29] stated the need of further approximation of ∂ n y to a continuous form on the triangulation as ∂ h n y h which may need extra calculation for the approximation of OCP. The error appearing in this interpolation will be contributed in the global error estimation of problem.
Instead of the standard FEM, we use the WEB-splines of degree n ≥ 2. WEB-splines are stable and their approximation order for solving elliptic PDE is O(h d ), where d is the degree of B-spline.
Furthermore, normal derivative of WEB-splines basis function is continuous and therefore the direct implementation of the normal derivative of state variable y h is feasible. In this way, an additional step can be avoided in the approximation of OCP.

Weighted Extended B-Splines
The representation and topological form of a simulation domain Ω may play a key role in choosing the best suited tool for related problems. Intensive research over the last three decades has proved B-splines are systematic techniques in approximation, computer aided design and graphics. B-splines as finite element basis can offer better solutions for both geometry description and numerical simulation [16]. B-splines are the piecewise continuous functions, usually represented by sequence of knots. Let (x τ , x τ+d+1 ) be a knot sequence, uniformly partitioning Linear, quadratic, cubic and higher order B-splines can be defined on this pattern. Stable and significantly more accurate approximated solutions are observed by using B-splines, especially in one dimension, as compared to classical finite element techniques. While for higher dimensions, (s ≥ 2) due to uniform support of B-splines, some limitations may arise when use them as finite element basis. Mainly, there are two reasons. Firstly, the use of Dirichlet boundary conditions is inconvenient as an approximation should be zero on the boundary and outside of a domain. The coefficients of the B-splines, whose supports have non empty intersection with boundary, must be zero to attain this task. This destroys the approximation order. Secondly, if some B-splines have very small intersection with a domain, the basis will not be uniform and become unstable near the boundary. This enlarges the condition number of system that results low convergence rate. These obscurities can be overcome by describing a new type of modified tensor product B-splines called weighted extended B-splines. WEB method was innovated and originated by Höllig et al. [30][31][32][33][34][35]. It has some major advantages. WEB method can be implemented on regular grids and there is no need of mesh generation as it is needed for the standard FEM. Basis functions are defined directly on grid. The domain retains its shape and is not disturbed by domain discretization or meshing. Boundary limitations are fully satisfied by this method. This method diminishes the size of the Galerkin system which ensures the accuracy of the approximated solution, fast convergence of solution algorithm, low condition number and less computation cost. For the sake of simplicity, we define here u v to be analogous to writing that u ≥ cv, where u and v are the functions and c is any constant independent of h.
More concerned are those B-splines whose support intersects a domain. These B-splines form the spline space B. Such B-splines are referred to as relevant B-splines. The set of indices of such B-splines is denoted as K. Additionally K is partitioned into inner and outer sets. The set of inner indices I is defined in such a way that the support of B-splines have at least one grid cell lying inside a domain.
This array corresponds to every outer index. The construction of extension coefficients, which couples the inner B-spline to the corresponding outer B-spline or vice versa, depends upon the value of Lagrangian polynomial [31]. The extension matrix [e i,o ] can be formulated as: where l λ is given in (16). This extension coefficient comprises the weighted extended B-splines defined as follows: x ∈ Ω as a positive weight function of order γ. The weight function ω should be continuous throughout the boundary and uniform on domain Ω, i.e., the weight function is sufficiently regular on a boundary. Ravashev et al. [36] studied R-functions in detail in order to construct weight functions.
However, for the case "max |B i (x)| 1", B i (x) cannot be used as a finite element basis. To normalize the basis, we have to multiply B i by balancing factor 1/ω(x i ), i.e., 1 The WEB space B web has less dimension than the dimension of the B-spline space B and holds all basic requirements of finite element basis functions.

Approximation of WEB Method
General results which are to be used in approximation of WEB-splines are stated. Höllig [30] gives the detailed descriptions and proofs of these results. The weight functions are considered to be k-regular and have order γ. The regularity of the continuous function plays a vital role in approximation results.

Some Useful Results
Quasi Interpolation [30]: In case of WEB-splines, the standard projector is defined by quasi interpolation as where Ξ i is the dual function for WEB-splines, also containing a finite number of cells in their support and is bounded uniformly. Inside a grid unit M, the norm of the standard projector is bounded, where M is the union of supports of all the relevant WEB-splines that intersect M ∩ Ω. The parameter c 1 is the constant that depends on domain, weight function and degree of the B-splines.
Bramble-Hilbert Estimate [30]: In scaled Ω, for orthogonal projection Q q upon the polynomial g with degree of at most "n" satisfies for 0 ≤ ξ ≤ θ ≤ n + 1 and c 4 is based on the domain and degree of B-splines.

Numerical Examples
The rational bezier curves are used to design the boundaries which are better to design and control the shape of the domain. Weight functions are the usual R-functions. Cubic B-splines are used to approximate the OCP. The solution behaves differently and has a variety of proportions.
The domain of a unit circle is taken with Ω ad = Ω. The function f = 1 − u, the matrix A = I and the parameter β = 1 is used. The solution and error plot of state equation are given in Figure 3. Figure 3a shows that the solution of the state equation is uniform, smooth and fully satisfies the Dirichlet boundary conditions, while the Figure 3b shows the error graph. The error is rigorous along the boundary. This behavior of error is due to the source function. Figure 4a shows the solution of the adjoint state equation. The solution does not show a good behavior along the boundary due to the normal vector ∂ n y. However the solution satisfies the boundary conditions of adjoint state equation. Figure 4b shows the error plot of adjoint state equation. The error plot shows distortion along the boundary while it is uniform inside the domain.    The solution of the state equation is smooth and satisfies the Dirichlet boundary conditions as shown in Figure 5a, while the error plot in Figure 5b shows that the solution has singularity along the corners.
The solution of adjoint state equation is not good enough along the boundary, shown in Figure 6a. Again the normal vector ∂ n y is responsible for this unhealthy behavior. However, the boundary conditions are satisfied. Relative errors are shown in Figure 6b. The error plot shows irregular behavior along the boundary while it is uniform inside the domain. Table 2 shows the error norms for u, state equation and adjoint state equation.  In WEB method, convergence depends upon the degree of the polynomial. Therefore, in Tables 1 and 2, state equation gives better approximation as compared to classical FEM, while due to the regularity issues, the adjoint state equation gives a similar approximation as classical FEM.

Conclusions
WEB-spline-based finite element method is used to approximate the OCP with outward normal derivative boundary observations. The cubic B-splines are used for the extension of B-splines because of their continuous derivatives and thus the good approximations are achieved. The WEB method is defined on a usual grid and proved itself a reduced computation method when there is a question mark on the complications of mesh generation. The error results for simulating OCP on two different domains show that the WEB method is converging. The error graph shows that the solution is not smooth along the boundary. However, by using hierarchical structure, this boundary problem can be resolved.